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I .  INTRODUCTION 


This  report  discusses  research  carried  out  by  Berke¬ 
ley  Research  Associates,  Inc.  under  contract  #N00014-83-C- 


i 

2157  with  the  Plasma  Fjhysics  Division  of  the  Naval  Research 
Laboratory.  The  repor^covers  the  performance  period  1  May 
1983  to  30  April  1984.  During  this  period,  an  extensive 
investigation  of  electron  beam  stability  was  carried  out  in 
the  context  of  an  easily  formulated  2-D  cold  electron  layer 
model.  This  work  was  carried  on  in  close  association  with 
NRL  personnel  in  conjunction  with  the  special  focus  program, 

*  Advanced  Accelerators.- 

^The  modified  betatron  has  been  the  subject  of  an  in¬ 
tense  investigation  by  NRL  staff  as  a  potential  high  current 
electron  accelerator.  Significant  progress  has  been  made  in 
the  understanding  of  the  operation  of  this  device  and  a  large 
experimental  effort  designed  to  test  this  understanding  is 
now  well  underway.  Much  of  the  work  to  date  has  been  direc¬ 
ted  toward  a  formulation  of  the  in  situ  beam  dynamics,  in¬ 


volving  equilibrium 


and  stability 


:  i  4 , 5  r6 


considerations , 


while  other  research  has  been  directed  toward  solving  the 

^ — 8 — c 

difficult  injection  problem  '  . 

C  Once  an  equilibrium  is  established,  which  in  general 


requires  constraints  on  beam  density  and  on  the  betatron 
field  index,  stability  questions  of  two  general  types  arise.  ^ 
First  is  the  question  of  the  stability  of  the  betatron 


oscillations  in  the  presence  of  small  field  errors,  i.e. 
the  problem  of  orbital  resonances.  These  resonances  are 
basically  single  particle  effects  caused  by  a  repeated  in- 
phase  interaction  of  a  particle  with  a  field  imperfection. 
However,  even  though  they  are  basically  single  particle 
effects,  resonances  can  set  a  limit  to  the  beam  current  in 
an  accelerator  since  the  betatron  frequency  of  one  particle 
is  affected  by  the  collective  fields  of  the  other  particles 
in  the  beam.  These  resonant  instabilities  were  studied 
under  an  earlier  contract  with  NRL,  N00014-81-C-2371 . 

The  second  broad  class  of  instabilities  occurs  in  a 
device  with  "perfect"  (azimuthally  symmetric)  applied  fields. 
These  generally  come  under  the  heading  of  collective  effects 
and  include  the  various  mechanisms  by  which  the  beam  may 
bunch,  kink,  or  simply  spiral  toward  the  chamber  wall. 

They  include  the  negative-mass-kink  mode  and  the  various 
resistive  wall  instabilities,  all  of  which  place  certain 
constraints  on  beam  current. 

The  various  collective  instabilities  in  the  modified 

betatron  were  first  studied  by  Sprangle  and  Vomvoridis5,  who 

characterized  the  different  modes  and  analyzed  the  dependence 

of  the  growth  rates  on  system  parameters.  Building  on  this 

9 

work,  Sprangle  and  Chernin  investigated,  under  somewhat 
more  general  conditions,  the  behavior  of  the  growth  rates  as 
current  was  raised,  keeping  the  beam  energy  spread  fixed. 


. 
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A  rather  unusual  and  unexpected  result  was  found,  namely, 
a  double-valuedness  in  the  current  versus  energy  spread  sta¬ 
bility  curve  and  the  appearance  of  two  disjoint  stability 
regions.  This  double-valuedness  was  attributed  to  the  ef¬ 
fects  of  self  fields  and  was  found  to  be  the  result  of  a 
competition  between  the  growth  and  stabilization  mechanisms. 
Under  some  conditions,  the  double  valuedness  could  even  be 
present  for  a  conventional  betatron  configuration. 

Prompted  by  this  novel  finding  and  mindful  that  the 
three-dimensional  nature  of  the  problem  made  inclusion  of 
self  field  effects  necessarily  approximate,  a  rigorous  study 
was  begun  of  a  two-dimensional  model,  including  correctly 
the  crucial  radial  and  azimuthal  particle  dynamics,  but 
neglecting  the  inessential  dynamics  in  the  axial  direction. 
This  model  has  the  advantage  that  all  self  field  effects 
could  be  included  exactly  in  the  formulation  of  the  eigen¬ 
value  problem  governing  the  small  signal  behavior.  This 
study  quickly  developed  in  several  interesting  directions, 
having  now  produced  potentially  useful  results  for  beam  be¬ 
havior  in  accelerators,  storage  rings,  and  microwave  tubes. 
It  is  the  results  of  this  study  which  are  reported  here. 

The  following  section  describes  the  model  we  have 
used,  discusses  its  essential  features,  and  summarizes  our 
findings.  Detailed  analysis  is  given  in  papers  attached  as 
Appendices.  Also  included  as  an  Appendix  is  documentation 


II.  TWO-DIMENSIONAL  LAMINAR  ELECTRON  LAYER  MODEL 


The  negative  mass  instability,  diocotron  instability, 
and  electron  cyclotron  maser  instability  all  depend  in  some 
way  on  the  existence  of  shear  in  the  fluid  velocity  in  the 
underlying  equilibrium  state.  Self  fields  in  an  intense 
beam  may  affect  the  magnitude  and  even  the  sign  of  the  shear 
in  a  beam;  hence  it  is  important  to  include  self  field  ef¬ 
fects  accurately  in  any  analysis  of  these  beam  instabilities. 

The  cold  fluid  equilibria  in  a  torus  are  quite  diffi¬ 
cult  to  calculate^.  The  stability  analysis  is  even  more  in¬ 
volved.  To  make  progress,  certain  simplifying  assumptions 
about  the  functional  form  of  the  self  fields  are  often  made  . 
These  assumptions  may  affect  the  dependence  of  the  shear  on 
the  beam  density.  One  is  therefore  motivated  to  consider  a 
"beam"  model  which  retains  the  important  effects  of  curvature 
while  allowing  an  exact  treatment  of  self  field  effects. 

Such  a  model — the  E-layer  or  Astron  model — is  illustrated  in 
Figure  1.  Here  a  laminar  layer  of  electrons  is  supported  by 
a  combination  of  radial  electric  and  axial  magnetic  fields 
between  two  coaxial  cylinders.  Here  the  cylinders  will  be 
taken  to  be  perfect  conductors,  though  resistive  wall  and 
other  boundary  effects  may  be  included  in  a  straightforward 
way,  as  discussed  in  Appendix  B. 

The  E-layer  model  has  been  the  subject  of  numerous 
studies^  Surprisingly,  however,  the  exact  linear 


stability  analysis  for  a  finite  thickness,  relativistic 
layer  does  not  seem  to  have  been  carried  out  previously. 

’so,  perhaps  as  a  consequence,  some  very  interesting,  poten 
tially  useful  features  of  the  model  appear  to  have  been  mis¬ 
sed  and  the  unifying  relationships^  among  the  various  insta 
bilities  seem  not  to  be  widely  appreciated. 

In  the  basic  equilibrium  state  illustrated  in  Figure 
1,  only  a  radial  electric  field  (Er)  and  an  axial  magnetic 
field  (Bz)  exist.  These  fields  are  in  general  due  to  a  com¬ 
bination  of  self  and  applied  fields.  The  equilibrium  is 
governed  by  the  following  set  of  three  equations  (MKS  units) 


y  v2 

-  —  °  =  -  [E  +  v  B  ] 
r  m  1  o  o  oJ 


r  ' 


-eVEo 


dB_ 
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These  equations  must  be  solved  subject  to  a  specified 
potential  difference  between  the  inner  and  outer  walls  and 
specified  total  magnetic  flux.  A  method  for  doing  this  and 
some  features  of  the  solution  are  discussed  in  Section  II  of 
Appendix  B.  Here  we  note  only  that  solutions  may  be  gener¬ 
ally  grouped  into  three  types,  depending  on  which  two  of  the 
three  terms  in  the  equation  of  force  balance,  Eq.  (1) ,  are 
most  important.  The  three  cases  are:  (I)  magnetic  force 


balances  centrifugal  force  (e.g.  Astron,  gyrotron,  particle 
accelerators),  (u)  electric  force  balances  magnetic  force 
(e.g.  crossed  field  microwave  devices) ,  and  (III)  electric 
force  balances  centrifugal  force  (e.g.  orbitron^4,  helio- 
tron^^) .  Our  analysis,  therefore,  has  wide  applicability 
and  many  standard  results  may  be  recovered  in  appropriate 
limits,  as  discussed  in  Appendix  B. 

One  of  the  most  important  equilibrium  quantities  is 
the  dependence  of  the  rotation  frequency  =  vQ/r  on  the 
energy  E  =  mc2y  -  e<p.  This  quantity  may  be  shown  to  be  given 


^o  =  *0  1~  Bo+2hH/Yo~l 

E  “me *Y0B£  [I+yJH^hTyJJ 


where  £  =  “p/^o'  wp  *s  tlie  Plasma  frequency,  and 


"  “  mc*B2Y3  ' 

o '  o 

When  is  a  decreasing  function  of  energy,  one  ex¬ 
pects  the  usual  negative  mass  instability  mechanism  to  oper¬ 
ate.  However,  we  note  the  possibility  of  reversing  the  sign 
of  d«0/dE  by  judicious  choice  of  applied  electric  field.  In 
fact,  as  detailed  analysis  shows,  the  negative  mass  insta¬ 
bility  may  be  completely  stabilized  by  choosing  +  2h  <  0. 
A  discussion  of  this  result  is  given  in  Appendix  A  and  a 
detailed  derivation  in  Appendix  B.  Another  noteworthy  fea¬ 
ture  of  Eq.  (4)  is  that  it  possesses  a  (negative)  extremum 


as  a  function  of  h,  at  h  =  y~  .  Since  hy£  is  the  ratio  of 
electric  to  centrifugal  forces  we  identify  those  devices  of 
Type  III  as  the  most  unstable  type  of  microwave  generators. 
This  result  suggests  that  the  orbitron  and  heliotron  might 
possess  some  advantages  over  other  microwave  tubes.  This 
possibility  is  presently  under  investigation. 

The  linear  stability  of  equilibria  described  by  Eqs. 
(1-3)  to  perturbations  in  the  r  and  6  directions  has  been 
studied  in  some  detail.  The  eigenvalue  problem  for  the  TE 
modes  has  the  form 

ra7(rSaf)  +  «  -  0  1 

where  4>  =  rEQ1  and  where  A  and  C  are  complicated  functions 
or  r,  the  azimuthal  mode  number  l,  and  the  eigenvalue  w. 

The  full  equation  is  derived  in  Appendix  B  where  the  features 
of  the  equation  are  also  discussed.  We  emphasize  that  Eq.  (6) 
is  fully  relativistic,  fully  electromagnetic,  and  contains 
all  effects  of  self  fields. 

Eq.  (6)  describes  two  basic  unstable  modes,  the  longi¬ 
tudinal  mode,  satisfying  w  -  luQ=0  and  the  transverse  mode, 
for  which  w  -  lwQ  =  ±  u>0(l  +  Y^h2)  2.  The  longitudinal  mode 
involves  significant  azimuthal  bunching  while  the  transverse 
mode  involves  mostly  radial  motion  of  the  layer.  In  the 
planar  limit,  it  may  be  shown  that  the  transverse  mode  re¬ 
duces  to  the  so-called  Buneman  instability^,  first  studied 
in  connection  with  the  early  magnetrons.  In  cylindrical 


geometry  the  transverse  mode  has  been  invoked  to  explain  the 

^  ^  ^  17 

operation  of  the  peniotron 

Our  studies  have  been  concentrated  on  the  longitudinal 
mode  which  includes  the  negative  mass,  cyclotron  maser,  and 
diocotron  instabilities.  For  a  thin  layer  it  is  possible  to 
obtain  an  analytical  approximation  to  the  growth  rate  for  the 
longitudinal  mode,  as  shown  in  Appendix  B: 


(Im(io)  ]  2  = 


(B2q+ 2h) 


(b++b_)  (1+Y^h2 ) 


+  OUi/R) 


where  t  =  layer  thickness,  R  =  layer  radius,  b+  =  wave  admit¬ 
tances.  (See  Appendix  B,  where  the  second  order  term  is  also 
given  explicitly.)  This  expression  for  the  growth  rate  has 
several  important  features:  The  growth  rate  (1)  scales  as 
(current)  ;  (2)  vanishes  when  +  2h  <  0;  (3)  has  a  maximum 
as  a  function  of  h  when  hy£  =  1;  and  (4)  vanishes  in  the 
planar  limit.  Features  (2)  -  (4)  are  novel  findings.  They 
are  expected  on  physical  grounds,  yet  they  have  not  been  pre¬ 
viously  discussed  in  the  literature.  Item  (2)  may  have  appli- 

18 

cation  to  low  energy  (~MeV)  storage  rings  .  Note  that  self 
field  effects  were  included  in  the  derivation  of  Eq.  (7) . 

The  stability  condition,  (2) ,  therefore,  is  insensitive  to 
beam  current,  azimuthal  mode  number,  or  container  geometry. 

Extensive  numerical  tests  of  Eq.  (7)  have  been  carried 
out,  using  a  computer  code  developed  to  solve  Eq.  (6)  through 
using  a  "shooting  method." 


In  the  tests  all  parameters  but  one 


were  held  fixed  and  the  dependence  of  the  growth  rate  on  the 
single  variable  parameter  was  determined  and  compared  with 
Eq.  (7).  The  parameters  tested  were  h,  v/y  »  t/R,  b/R,  a/R, 

Y  ,  and  l.  Good  agreement  was  found  in  all  cases;  details 
are  given  in  Appendix  B.  The  computer  code  used  to  solve 
the  eigenvalue  problem  is  documented  in  Appendix  C. 

In  addition  to  these  tests,  a  search  was  carried  out 
for  evidence  of  the  double  valuedness  in  the  current  vs. 
energy  spread  curve,  the  prediction  of  which  was  the  origi¬ 
nal  motivation  for  this  work.  Despite  much  effort,  no  such 
double  valuedness  was  uncovered.  The  search  is  complicated, 
however,  because  it  is  difficult  to  develop  a  set  of  param¬ 
eters  for  layers  with  the  same  energy  spread,  but  different 
currents  since,  for  large  currents,  the  energy  spread  is  a 
function  of  the  current.  In  fact,  current,  energy  spread, 
and  beam  thickness  are  connected  in  a  way  that  is  absent  in 
the  equilibrium  used  in  reference  9.  The  reason  is  that  in 
reference  9,  betatron  oscillations  are  assumed  to  "fill  in" 
the  beam  to  the  size  and  shape  assumed  when  calculating  the 
self  fields.  Inclusion  of  betatron  oscillations  in  the 
E-layer  model  may  help  to  carry  out  an  effective  comparison 
more  readily,  but  this  is  not  easily  done.  A  definitive 
resolution  of  the  question  on  the  existence  of  double-valued- 
ness  must  await  further  analysis,  though  numerical  particle 
simulations  may  also  be  helpful  in  shedding  some  light  on 
this  question. 


III.  SUMMARY  AND  RECOMMENDATIONS 


A  detailed  study  of  a  two-dimensional  electron  layer 

confined  between  coaxial  cylinders  has  been  carried  out  with 

the  intention  of  understanding  the  radial  and  azimuthal  beam 

dynamics  in  high  current  accelerators.  Though  no  confirma- 

g 

tion  could  be  found  for  the  predicted  double-valuedness  in 
the  current  versus  energy  spread  curve,  several  interesting 
and  potentially  useful  features  of  electron  beam  behavior 
in  curved  geometries  were  discovered.  Specifically,  it  was 
found  that  by  appropriate  choices  of  applied  electric  field, 
the  growth  rate  of  the  longitudinal  mode  could  either  be 
maximized,  as  might  be  desirable  in  certain  microwave  devi¬ 
ces,  or  eliminated  altogether,  as  would  be  desirable  in 
electron  accelerators  and  storage  rings.  In  addition,  a 
long-standing  puzzle  concerning  the  survival  of  a  finite 
growth  rate  for  the  negative  mass  instability  in  the  planar 
limit  was  resolved. 

There  are  several  areas  where  further  work  would  be 
useful.  First,  additional  research  is  needed  to  resolve  the 
question  of  double  valuedness  discussed  above.  The  problem 
is  a  difficult  one,  however,  and  computer  simulations  may 
be  needed  to  resolve  this  complex  issue.  Second,  the 
effects  on  beam  behavior  of  other  applied  magnetic  or  elec¬ 
tric  fields  is  also  an  important  area  for  further  investiga¬ 
tion.  In  particular,  inclusion  of  a  toroidal  field  in  the 


E-layer  model  is  straightforward  and  leads  to  a  fourth  order 
eigenvalue  problem  whose  features  are  important  to  study. 
Careful  examination  of  this  problem  may  lead  to  a  more  thor¬ 
ough  understanding  of  the  effect  of  the  toroidal  field  on 
beam  stability  in  the  modified  betatron.  And  finally,  we 
mention  that  the  effect  of  strong  focusing  fields  on  these 
instabilities  is  an  important  topic  for  future  study,  since 
it  has  been  shown  that  use  of  such  fields  may  have  certain 
beneficial  effects  on  the  beam  behavior.  Even  these  fields 
may  be  included  in  the  E-layer  model  if  0-dependence  is 
allowed  in  the  equilibrium  state. 
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li  is  shown  (hat  the  negative  mass  instability  in  a  rotating  relativistic  electron  layer  may  be 
stabilized  by  a  radial  dc  electric  field  of  a  suitable  magnitude  The  stabilization  mechanism  is 
independent  of  the  beam  velocity  spread,  and  is  insensitive  to  the  beam  current,  the  con¬ 
tainer  geometry,  or  the  azimuthal  mode  number.  A  simple  stability  criterion  is  given. 

PACS  numbers  52  60  4  h.  29  20  -c.  47.75  4  f,  52.3S.P> 


The  negative  mass  instability1  poses  a  major  ob¬ 
stacle  to  the  development  of  high-current  cyclic  ac¬ 
celerators.  Various  methods  of  stabilization  have 
been  proposed  and  analyzed.  Notably,  the  effects 
of  a  moderate  beam  angular  velocity  spread  and  be¬ 
tatron  oscillations  have  been  considered  2-6  For  a 
betatron,  the  addition  of  a  toroidal  magnetic  field7 
has  been  shown  to  reduce  the  instability  growth  rate 
considerably51  and  for  the  Astron,  the  proximity  of 
the  container  walls  to  the  relativistic  electron  layer 
(Flayer)  stabilizes  the  lower  azimuthal  modes.1 

Jr  this  Letter,  we  show  that  by  imposing  a  nega¬ 
tively  biased  radial  electric  field  of  a  suitable 
strength,  the  negative  mass  instability  may  be 
suppressed.  This  stabilization  differs  from  all  previ¬ 
ously  known  mechanisms  in  that  it  is  effective  even 
for  a  very  cold  beam,  it  does  not  require  a  toroidal 
magnetic  field,  nor  is  it  sensitive  to  the  container 
geometry,  the  beam  current,  or  the  toroidal  mode 
number.  The  simple  stability  criterion,  given  in  Eq. 
(10)  below,  does  not  seem  to  be  very  stringent  for 
electron  beams  in  the  megaelectron  volt  range. 

Our  finding  is  based  on  an  analytic  treatment  of 
the  stability  of  the  £  layer  situated  in  a  configura¬ 
tion  similar  to  the  Astron,  which  has  been  shown* 
to  include  all  essential  features  of  the  negative  mass 
instability.  We  limit  our  study  to  a  highly  ordered 
beam  whose  unperturbed  orbits  are  concentric  cir¬ 
cles.  Such  a  beam  should  yield  the  most  pessimistic 
prediction  as  far  as  the  beam  stability  is  concerned; 
hence  our  analysis  is  conservative.  The  simplicity 
of  the  assumed  equilibrium  orbits  allows  the  linear 
stability  theory  to  be  formulated  exactly,  including 
alt  ac  and  dc  space-charge  effects,  all  relativistic  ef¬ 
fects,  and  all  electromagnetic  effects,  for  general 
equilibrium  profiles.  As  we  shall  see,  our  disper¬ 
sion  relation  reproduces  the  standard  results  in  the 
appropriate  limits.  For  example,  the  diocotron  in¬ 
stability  is  recovered,  and  the  negative  mass  insta¬ 
bility  removed,  in  the  planar,  nonrelativistic  limit. 


Consider  a  cylindrical  £  layer  with  radial  density 
profile  no(r)  which,  in  equilibrium,  circulates  con¬ 
centrically  with  azimuthal  velocity  v0(r)  -  Ov0(r) 
m0ruo(r)  under  the  combined  action  of  an  axial 
magnetic  field  B0  ~iB0(r)  and  radial  electric  field 
Eo*r£o(r).  These  fields  include  both  the  self- 
fields  and  the  externally  imposed  fields.  We  as¬ 
sume  that  the  £  layer  is  located  between  two 
cylindrical  conductors  of  inner  and  outer  radii  a  and 
b.  respectively,  and  that  there  is  no  axial  motion 
nor  axial  variation  in  either  the  unperturbed  or  the 
perturbed  states 

The  governing  equations  for  the  equilibrium  read 


y0vl/r-  -  (e/m0)(£0  + voR0). 

(1) 

dBcJdr  *  — 

(2) 

r_1  d(rE0)/dr «  n0e/f0 

(3) 

Here,  e  and  m0  are  respectively  the  electron  charge 
and  rest  mass,  po  and  <o  are  the  free-space  permea¬ 
bility  and  permittivity,  and  y0“  d  -  vl/c7)~m  is 
the  relativistic  mass  factor  with  c  being  the  speed  of 
light.  Once  the  electron  density  n0(r),  the  total 
electrostatic  potential  difference  between  r  «  a  and 
r  -  b,  and  the  magnetic  flux  are  specified,  the  un¬ 
perturbed  fields  v0(r).  £(>(').  and  fi„(r)  are  to  be 
solved  from  (1)  to  (3)  to  yield  a  self-consistent 
equilibrium  solution. 

We  next  consider  a  small-signal  perturbation  on 
such  an  equilibrium.  All  perturbations  are  assumed 
to  vary  as  fir)  exp(  iut  -  HO),  where  /  is  the  azimu¬ 
thal  mode  number  and  u>  is  the  (complex)  eigenfre- 
quency  to  be  determined.  In  the  absence  of  axial 
variation  and  of  axial  motion,  the  TM  modes  and 
the  'jTE  modes  are  decoupled.  The  nontrivial  com¬ 
ponents  of  the  rf  electromagnetic  fields  are  £,,  £,. 
and  B ,  for  the  TE  modes.  The  Maxwell  equations, 
the  Lorentz  force  law,  and  the  continuity  equation 
may  then  be  linearized  and  combined  to  yield  the 
following  second-order  ordinary  differential  equa- 
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lion  for  4> »  rE, 

In  this  equation 

ft  ■<!»- lu>o(r),  Sp  m uif/wlD,  p»  1  -  l2c 2/u>2r2  +  Vpft  2/u>2, 

eJno/yo^o*0'  D  m  PQ  —  ft5/«*»o,  P  “  yo^  1  +  6 ).  0  *  6  +  v^/wq, 
i)d“ SD(uo/u)P/yl.  A  »  -tj DU/pr})(  \  +  Npfild/lui0), 

/30“  vq/c,  0  »  (u2/c2r)(  1  -  rje/p  +  Nphl)- 

In  the  definition  of  a  prime  denotes  a  derivative  with  respect  to  r.  and  h  is  proportional  to  the  equilibrium 
electric  Held  and  is  defined  by 

h  ■  -  erE(Jmnylvl.  (5) 

Note  that  h  is  positive  (negative)  if  the  equilibrium  electric  field  points  radially  outward  (inward).  The 
eigenvalue  <•>  is  determined  by  solving  (4)  subject  to  the  boundary  conditions  d>  *  0  at  r  «  a  and  at  r  -  b. 

Equation  (4)  is  completely  general  and  of  wide  applicability.  It  governs  the  small-signal  stability  properties 
of  various  devices  including  the  Astron,J,4r9  gyrotron,10  orbitron,11  and  cross-field  microwave  devices,12 
depending  on  the  parameters  of  the  electron  beam  as  long  as  the  equilibrium  states  are  modeled  by  Eqs. 
(l)-(3).  A  detailed  comparative  stability  study  of  various  types  of  equilibrium  will  be  given  elsewhere.  For 
the  present  purpose,  we  restrict  ourselves  to  an  £  layer  with  uniform  density  n0  extending  from  r«r,  to 
r-r2.  The  £-layer  thickness  r«r2-ri  is  assumed  to  be  much  less  than  the  mean  radius  R  We  shall  use 
t/R  as  an  expansion  parameter.  Furthermore,  we  assume  that  |ft|  «  w„.  a  condition  readily  satisfied  by 
the  negative  mass  mode.'"10 

The  instability  growth  rate  w,  may  be  analytically  derived  from  Eq.  (4)  for  a  thin  £  layer  by  expanding 
about  the  singularity  ft  -  0  in  the  complex  r  plane.  To  two  orders  in  t/R,  it  is  given  by 


A+ +6- 


/t  (Po  +  2h )  fir*  2 

R  (1  +ylh2)  R2 


where  6+  ib.  )  is  the  normalized  wave  admittance  at  the  outer  (inner)  edge  of  the  £  layer.13  and 

.  I _ ! _ 

4  y$(l  +  6  )2(1  +ylh2) 

*\h  100+26  +  (1  +  6  )21  +  2  ^  yS(/3g+26)2-yg(l+y$62)(|SS+26)2 

(l«6  "0  J 

1  ( 1  +  h  )wf/<i>l  |  b  +  -  b  _  |1? 

~  2  yUl  +  yfr2)  I  6+  +6-  |j 

The  first-order  term  (in  t/R  )  of  (6)  describes  | - — 

the  negative  mass  effect1  and  dc  field  effects,  while  (b)  A  more  stringent  test  on  the  validity  of  (6)  is 
the  second-order  term  includes  the  diocotron  ef-  to  consider  the  planar  geometry  limit.  In  this  limit. 

feet14  and  finite-thickness  stabilization.3  The  we  let  R  —  «,  /—  »,  <u0  —  0,  but  require  that  v0. 

derivation  of  (6)  will  be  given  elsewhere.  Its  validi-  ky  "l/R,  £0,  and  t  remain  finite.  Then  h  —  <» 

ty  may  be  tested  as  follows:  from  (5)  and  the  first  term  of  the  right-hand  side  of 

(a)  If  the  beam  is  infinitesimally  thin,  and  if  we  (6)  tends  to  zero,  consistent  with  the  notion  that 

ignore  the  dc  electric  field  by  setting  h  -0.  we  then  there  is  no  negative  mass  instability  in  a  planar 

recover  from  (6)  the  well-known  dispersion  rela-  geometry.  Using  only  the  last  term  of  (6).  we  then 

tionship  *»}  a  (/t/R  )w}&l/(b+  +  b  _ )  for  the  neg-  obtain 

ative  mass  instability  for  the  Astron  geome-  ,  ,,  ,  .  .  ,  ,  ,,  .  .  , 

,rwM.«.l0  =  (*,V/4y&62)o,;M-  U,V/4y5)w;/u,2 


15-2 


*  A/.  A  • . 
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which  agrees  with  the  well-known  growth  rate  for 
the  diocotron  instability  for  a  sheei  beam  14  In  writ¬ 
ing  the  last  expression,  we  have  used  the  self- 
consistent  equilibrium  condition  £o+v0*f?o“0 
(for  a  planar  sheet  beam)  in  (5)  and  defined 

*  klRt/rw0>o  This  agreement  with  previously- 
known  results  adds  to  our  confidence  in  the  disper¬ 
sion  relationship  (6).  especially  with  regard  to  the 
effects  of  dc  self-fields  Recent  work  on  the 
diocotron  instability  is  reported  by  Tsang  and 
Davidson  1J 

For  a  thin  £  layer  with  sufficiently  high  energy 
( ;>  1  MeV),  the  last  term  of  (6)  may  be  ignored. 
The  dispersion  relationship  may  then  be  approxi¬ 
mated  by 


+6. 


/r  |  (06+26) 

R  (1  +  y^?) 


Thus,  the  sufficient  condition  for  stability  is 


h  <  -pl/2 


(9) 


for  the  usual  case15  b*+b_  >0.  This  stability- 
condition  (9),  together  with  the  definition  of  h  in 
(5),  implies  that  a  sufficiently  strong,  radially  in¬ 
ward  electric  field  may  render  the  relativistic  elec¬ 
tron  beam  stable  against  the  negative  mass  instabili¬ 
ty.  This  stabilization  is  independent  of  the  beam 
velocity  spread  or  betatron  motion,  and  since  its 
derivation  has  already  taken  into  consideration  the 
dc  self-field  effects,  the  criterion  (9)  is  not  restrict¬ 
ed  to  a  low -current  beam  16  In  the  event  that  the 
externally  imposed  electric  field  exceeds  the  self- 
fields.  the  stability  condition  (9)  may  be  rewritten 
as 


|ed»l  >  ’T’^o>o  ln(b/o  ),  (10) 


where  |<*d>|  is  the  externally  imposed  potential 
difference  (in  kiloelectronvolts)  between  r  «■  a  (the 
cathode)  and  r  -  b  (the  anode). 

As  an  example,  take  R  *  100  cm,  b  -  c  *=  4  cm. 
Then,  according  to  (10),  a  I -MeV  electron  beam 
would  be  stable  against  the  negative  mass  instability 
if  the  inner  conductor  is  negatively  biased  at  a  volt¬ 
age  greater  than  200  keV  with  respect  to  the  outer 
conductor. 

A  partial  explanation  of  the  stability  condition  (9) 
may  be  given  in  terms  of  the  single  particle  motion 
in  an  externally  imposed  field  £0  and  Bb  Let  «  be 
the  total  energy  (kinetic  and  potential)  of  an  elec¬ 
tron.  One  may  easily  deduce  from  (1)  that 
dwrjdt  m  (dtujdr  )dr/dta  (0j,+  2 h  )/(!-+■  y\h J)  if 
the  self-fields  of  the  electron  layer  are  neglected. 
Thus,  the  effective  azimuthal  inertia1-9  of  a  rotating 


relativistic  electron  may  be  converted  from  negative 
to  positive  if  h  is  less  than  It  should  be 

stressed,  however,  that  the  stability  condition  (9)  is 
derived  from  collective-mode  considerations  which 
include  both  the  ac  and  the  dc  self-fields 11 

In  summary,  this  Letter  presents  a  novel,  robust 
method  to  suppress  a  major  instability  in  circular  ac¬ 
celerators.  Technical  aspects  such  as  fabrication, 
beam  injection,  and  beam  retrieval  remain  to  be 
studied.  A  more  refined  analysis  may  be  needed  to 
examine  the  possible  occurrence  (if  any)  of  residual 
instabilities  The  stability  criterion  (9)  may  be  test¬ 
ed  on  several  currently  operating  devices. 
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•*If  the  container  wall  is  lossy,  the  wave  admittance 
(P*  +  6.)  would  be  complex.  The  electron  beam  may 
then  be  subject  to  resistive  wall  instabilities  even  if  the 
stability  criterion  (9)  is  satisfied.  The  resistive  growth 
rate  and  the  negative  mass  growth  rate  scale  differently, 
however.  We  wish  to  thank  A.  M  Sessler  (private  com¬ 
munication)  for  reminding  us  of  the  importance  of  resis¬ 
tive  wall  instabilities,  and  for  furnishing  an  argument 
supporting  our  conclusion  on  the  stabilization  mecha¬ 
nism. 


‘’it  is  of  some  interest  to  note  that  the  negative  mass 
factor  dwo/dt  is  maximized  with  respect  to  h  when 
h  -  \/yl  Thus,  according  to  Eqs  (1)  and  (5).  the  £ 
layer  is  most  unstable,  and  is  therefore  most  likely  to 
yield  radiation,  if  its  equilibrium  rotation  is  iohly  sup¬ 
ported  by  a  radially  outward  electric  field  Reference  11 
reported  a  potent  radiation  source  of  this  type  More¬ 
over.  since  the  dispersion  relationship  (8)  is  applicable 
for  arbitrary  combinations  of  £0and  Bo  and  for  arbitrary 
energy  of  the  electron  beam,  it  provides  a  ready  compar¬ 
ison  of  the  "potency"  among  various  microwave  devices 
such  as  the  gyrotron.  orbitron.  heliotron.  and  cross-field 
devices  (if  the  small-signal  growth  rate  is  used  as  a  cri¬ 
terion)  Further  discussions,  as  well  as  the  confirmation 
of  the  stability  criterion  (9)  by  a  numerical  integration  of 
(4).  will  be  reported  elsewhere 
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I.  Introduction 


There  exists  by  now  a  truly  enormous  literature  on  the  subjects  of  the 
equilibrium  and  stability  of  layers  of  charged  particles  in  various 
geometries.  The  earliest  studies1-2  were  conducted  in  order  to  understand  the 
“slipping  stream"  or  diocotron  effect  in  the  operation  of  the  first 
magnetrons.  Somewhat  later  the  importance  of  curvature  effects  was  realized 
when  beams  became  relativistic,  as  in  particle  accelerators;  the  resulting 

O  Q 

“negative  mass”  instability0’’  completely  dominates  the  planar  beam  diocotron 

effect  at  sufficiently  high  energies  (only  a  few  tens  of  keV  in  many  practical 

devices  of  interest).  The  negative  mass  effect  was  also  investigated  in 

connection  with  some  controlled  fusion  research  devices10-12  and  other 
13 

machines. 

In  recent  times  there  has  been  renewed  interest  in  high  power,  high 
efficiency  microwave  devices  as  well  as  in  accelerators  capable  of  high 
current  operation.  Spurred  by  the  discovery  of  the  electron  cyclotron  maser 
(or  gyrotron)  effect,14-16  research  in  the  field  of  short  wavelength,  high 
power  microwave  devices  requiring  no  slow  wave  structure  has  been  vigorously 
pursued.12  Operation  of  these  new  devices  depends  fundamentally  on  the 
negative  mass  effect  as  enhanced  by  a  synchronism  of  the  particles'  angular 
motion  with  the  temporal  and  angular  variation  of  a  “cold"  waveguide  mode. 

This  enhancement,  though  not  called  the  maser  effect,  was  first  noted  in  the 
classic  work  of  Briggs  and  Neil10  and  has  been  further  elaborated  in  ref.  18. 

The  acceleration  of  large  currents  of  electrons  is  a  formidable  problem 
which  has  also  received  considerable  attention  recently.  In  cyclic  devices  it 
is  possible  to  construct  high  current  equilibria  in  many  cases1^-21  but 
various  instabilities,  the  negative  mass  instability  prominent  among  them,  may 
limit  achievable  currents  to  smaller  values  than  equilibrium  considerations 


alone  would  suggest.  In  these  devices  the  effects  of  self  fields  on  the 
negative  mass  instability  and  on  the  stabilization  mechanism  become  important 
to  consider;  in  one  device  these  effects  have  been  predicted  to  lead  to  a 
peculiar  double-valued  feature  in  the  current  vs.  energy  spread  stability 
curve.  It  was  this  particular  result  which  initially  prompted  the  study  of 
self  field  effects  reported  here  for  a  simpler  (2D)  model. 

The  model  we  consider  consists  of  a  layer  of  charged  particles  (we  will 
think  of  them  as  electrons  but  ions  may  be  trivially  substituted)  moving  in 
circles  about  a  common  axis,  as  shown  in  Fig.  1.  This  restriction  to  laminar 
or  "cold"  electron  flow  will  tend  to  overestimate  actual  growth  rates  of  the 
modes  we  study  since  the  effects  of  betatron  oscillations  and  of  nonaxis- 
encircling  particles  are  stabilizing.  An  analysis  including  these  effects  is 
properly  done  in  phase  space  using  the  correct  equilibrium  distribution 
function;  such  an  analysis  is  significantly  more  complex  than  that  given  here 
and  may  be  carried  out  only  in  an  approximate  way.  The  laminar  flow  case 
includes  all  essential  physics  and  has  great  simplicity,  allowing  an  exact 
treatment  of  the  linearized  problem,  to  recommend  It.  The  equilibrium, 
discussed  in  the  following  section,  is  supported  by  a  combination  of  self  and 
externally  applied  radial  electric  and  axial  magnetic  fields.  We  Impose  no 
a  priori  restriction  on  the  relative  magnitudes  of  the  three  terms — 
centrifugal  force,  electric  force,  v  x  fc  force — in  the  equilibrium  force 
balance.  This  configuration  is  thus  a  reasonable  model  for  the  Astroo,23 
magnetron,  gyrotron,  orbitron,2*  peniotron25  and  heliotron26  and  includes 
correctly  the  crucial  radial  and  azimuthal  particle  dynamics  found  In  particle 
accelerators.  This  "E-layer"  model  has  the  virtue  that  the  linear  stability 
problem  may  be  formulated  exactly,  for  arbitrary  particle  energies,  including 
all  effects  of  self  fields  and  all  effects  of  relativity  and 


electromagnetism.  It  is  perhaps  a  bit  surprising  that,  despite  the 
venerability  of  the  topic,  this  exact  formulation  has  not  been  carried  out 
earlier. ^  The  eigenvalue  problem  governing  the  stability  of  the  E-layer  16 
derived  and  analyzed  in  Section  III,  below. 

The  desirability  of  a  completely  general  treatment  encompassing  many 
devices  and  a  large  parameter  space  is  related  to  the  ease  with  which  various 
familiar  results  for  special  cases  may  be  recovered  in  appropriate  limits.  A 
dispersion  relation  for  the  so-called  longitudinal  mode  obtained  analytically 
in  the  thin  beam  limit  from  the  ordinary  differential  equation  governing  the 
RF  field  reproduces  all  standard  special  case  dispersion  relations  (negative 
mass,  electron  cyclotron  maser,  diocotron)  in  a  straightforward  way;  in  doing 
so,  a  puzzle  is  resolved  concerning  the  survival  of  a  finite  negative  mass 

growth  rate  in  the  planar  limit  and  a  method  suggests  itself  on  how  either  to 

28 

maximize  or  eliminate  altogether  the  negative  mass  instability  growth. 
Extensive  testing  of  the  dispersion  relationship  against  a  numerical  solution 
of  the  eigenvalue  problem  shows  excellent  agreement,  as  reported  in  Section 
III.  The  classical  diocotron  dispersion  relation  is  discussed  and  recovered 
from  our  formalism  in  Section  IV. 

Our  treatment  of  the  transverse  mode,  which  has  been  invoked  in  the 

25 

theory  of  the  peniotron  device  is  somewhat  less  comprehensive.  This  mode  is 
unstable  only  when  electron  motion  is  synchronous  with  a  cavity  mode  and  so  is 
of  a  more  specialized  nature  than  the  negative  mass  or  diocotron  instabilities 
which  arise  from  intrinsic  properties  of  geometry  and  shear  flow,  not  from 
interaction  with  external  structures.  Still,  the  utility  of  the  transverse 
mode  for  microwave  generation  may  not  yet  have  been  fully  exploited.  We 
comment  on  some  features  of  this  mode  in  Section  V. 
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II.  Cold  E-Layer  Equilibrium 


We  consider  an  Idealized  model  of  an  electron  layer  in  which  all 
electrons  circle  a  common  axis  as  shown  in  Fig.  1.  We  will  neglect  the 
effects  of  betatron  oscillations  and  of  axial  motion  in  both  perturbed  and 
unperturbed  states.  Inclusion  of  betatron  oscillations  or  axial  perturbations 
is  expected  to  have  a  stabilizing  Influence  on  the  collective  modes  we  will 
study.  We  further  assume  that  in  general  the  layer  is  enclosed  in  a  coaxial 
waveguide  with  smooth,  perfectly  conducting  walls  at  r«a  and  r-b,  as  shown  in 
the  figure.  Absence  of  either  or  both  walls  results  only  In  a  change  of 
certain  boundary  conditions. 

Though  this  model  is  simple  to  describe  it  is  surprisingly  rich  in 

content.  Depending  on  beam  and  geometric  parameters  and  types  of  applied 

fields,  it  may  be  taken  as  a  good  description  of  the  Astron,23  gyrotron,17 
24  25  2ft 

orbitron,  peniotron,  heliotron*0  and  cross  field  microwave  devices  and  may 
also  be  of  some  Interest  in  accelerator®**  and  space  physics  stability 
problems  as  well  as  in  the  theory  of  magnetic  insulation. 

The  equations  governing  the  equilibrium  quantities  v  (r)0,  E  (r)r, 

A 

Bo(r)z  for  a  specified  density  profile  nQ(r)  are 


■  —  [E  +  v  B  I 
r  m  o  o  o* 


(1) 


I  i 

r  dr 


<rEo) 


-en  /e 
o  o 
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dB 

o 

-j —  *  u  en  v  , 
dr  Mo  o  o’ 


(3) 


where  -e  and  m  are  the  electron  charge  and  mass,  v 
2  -  1/2  0 
factor,  (l-(v  /c)  )  ,  c  is  the  speed  of  light  in 


is  the  usual  relativistic 

vacuum,  and  e  and  u  are 
o  o 


Che  permiccivicy  and  permeability  of  free  space.  In  addition  to  the  density 

profile,  ve  take  as  given  the  total  electrostatic  potential  difference  and  the 

total  magnetic  flux  contained  between  r*a  and  r*b. 

Using  these  conditions  the  solution  to  Eq.  (2)  is  immediate  and  Eqs.  (1) 

and  (3)  may  be  combined  to  yield  a  single  differential  equation  for 

u  =  B  Y  3  v  y  / c: 
o 1  o  o  o 


/ 

u 


- ~  y 

1  +  h 


2 

o 


h 

-] 


(A) 


2  2  2  2 

where  £  -  wp/wo,  o>  ■  e  n0/®Y0t0  is  the  plasma  frequency  of  the  layer, 
u>o  “  vo/r,  and  a  prime  denotes  d/dr.  In  Eq.  (4)  we  introduce  the  Important 
quantity  h,  defined  by. 


h 


erE 


2  2  3 
me  6  y 
o  'o 


C5) 


2 

which  is  (!/yo)  times  the  ratio  of  the  electric  force  to  the  centrifugal  force 
experienced  by  an  electron  in  equilibrium  at  radius  r. 

Equation  (4)  must  be  solved  subject  to  an  initial  condition  u(rj) 
satisfying  Eq.  (1)  which  we  rewrite  in  the  dimensionless  form 


u2  -  oe(u2+1)1/2  +  aBu 


(6) 


2 

where  Og  ”  erEo/mc  and  •  erBo/mc.  aE  is  known  at  r^  but  <*g  is  not.  The 

entire  equilibrium  problem  then  reduces  to  choosing  aB(r.)  such  that  the 
resulting  total  flux 


b 

♦o  “  2irJ  dr  rBo 

a 


(7) 
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is  a  specified  number.  A  numerical  method  for  doing  this  is  described  in 
Appendix  A.  Here  we  briefly  consider  the  character  of  some  of  the 
solutions.  The  equilibrium  thus  constructed  will  be  used  in  the  numerical 
solution  of  the  eigenvalue  problem  in  Section  III. 

r^),  Eq.  (6)  may  have  0,  1,  or  2  real  valued 
solutions  for  u(r^),  depending  on  the  values  of  otg  and  oB;  without  loss  of 
generality  we  take  Og(r^)  >  0.  The  situation  in  the  a^-  plane  is  then 
depicted  in  Fig.  2.  By  a  positive  root,  denoted  by  (+)  in  the  figure,  we  mean 
a  right-handed  rotation  about  $  (clockwise,  looking  in  the  direction  of  $). 

The  equation  for  the  boundary  curve  C  in  Fig.  2  is  the  condition  for  Eq.  (6) 
to  have  a  double  root: 


Given  OgCrj)  and  Og( 


where  u*  -  (8/3)1/2sinh[^  sinh"1  [^(3/2)3/2aB)  ] .  On  C,  1  +  h  -  0. 

There  are  three  basic  types  of  equilibria  found  in  devices  of  practical 
interest:  Type  I,  In  which  a  magnetic  field  is  used  to  balance  centrifugal 
force  (ag  is  small,  e.g.  Astron,  gyrotron,  particle  accelerators);  Type  II,  in 
which  electric  and  magnetic  forces  balance  (Inertia  is  small,  e.g.  crossed 
field  microwave  devices);  and  Type  III,  in  which  an  electric  field  is  used  to 
balance  centrifugal  force  (afi  is  small,  e.g.  orbitron,  heliotron).  Type  I  may 
have  only  "+"  roots  while  Types  II  and  III  may  have  either  or  roots. 

In  the  following  analysis,  however,  we  make  jw  assumption  regarding  the 
relative  magnitudes  of  the  electric  field,  magnetic  field,  and  inertia  terms. 

It  is  reasonable  to  ask  whether  stability  considerations  favor  one  type 


of  equilibrium  over  another  for  a  layer  of  given  kinetic  energy.  To  our 
knowledge  there  has  been  no  definitive  answer  to  this  fundamental  question. 


For  radiation  source  applications  one  arguably  sight  want  the  most  unstable 
configuration  while  for  other  applications,  like  accelerators,  one  would  want 
the  most  stable  one.  We  proceed  to  examine  the  stability  of  a  general 
equilibrium  descrlbable  by  Eqs.  (1)  -  (3)  in  order  to  address  this  question, 
our  goal  being  a  complete  parametric  study  of  the  small  signal  behavior. 

III.  Stability  of  the  Layer 

We  consider  perturbations  on  the  equilibrium  described  in  the  previous 
section.  The  equations  governing  the  layer  are,  simply 


(ft  +  “ 

-  +  Sat) 

ID 

(9) 

$•£  *  -en/e 

o 

(10a, b) 

\  l!  -  "oen? 

$•8  -  o 

(lla.b) 

where  n  is  the  number  density. 

Our  discussion  will  be  limited  to  consideration  of  perturbations  in 

at  ^ 

the  r  and  0  directions,  for  reasons  described  below.  Writing  all 
quantities  T  as 

T(r,©,t)  -  To(r)  +  e/^-J  ei(wt_£0)f t(£) (r ;u)  (12) 

where  ^(r)  is  the  equilibrium  value,  and  retaining  only  terms  linear 
ft) 

in  '  one  finds  that  the  linearized  versions  of  Eqs.  (9)  -  (11)  decouple 
into  two  sets  governing  {vlr,vl0,Elr .E^.B^.n^  and  "Mch 

we  Identify  as  TE  modes  and  TM  modes  respectively.  (Rere  and  below  we 


write  T1^(r;w)  simply  as  for  all  first  order  quantities.)  The  TM  modes 
■ay  easily  be  shown  to  be  neutrally  stable  for  perfect  conductor  boundary 
conditions  (damped,  for  resistive  wall  boundary  conditions).  They  represent 
simple  oscillations  of  the  electrons  along  the  z  direction  in  response  to  the 
cavity  mode  fields;  the  equilibrium  model  does  not  provide  any  free  energy  to 
excite  the  TM  modes.  The  TE  modes,  on  the  other  hand,  are  potentially 
unstable.  If  we  had  allowed  a  finite  axial  wavelength  for  the  perturbed 
quantities  the  TM  and  TE  modes  would  be  coupled,  but  it  has  been  found10  that 
the  effects  of  finite  axial  wavelength  are  stabilizing  for  the  coupled  case. 
Consequently  we  focus  attention  exclusively  on  the  TE  waves. 

The  Euler  equation  (9),  upon  linearization,  gives  the  fluid  response  in 
terms  of  the  RF  fields  of  the  TE  wave: 


'ln,tr~  ”0^10 


®Y, 


(E.  +  v  B,  1 
‘  lr  o  lzJ 


(13a) 


“oQvlr'1Bv10 


3  E10 


(13b) 
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where  we  have  defined  n  -  ,  P  *  y  (1-Hi),  Q  -  v7  /w  +  h.  and 

o  o  00* 

▼  '=  37  vq.  Using  the  solution  of  Eqs.  (13a, b)  in  the  definitions  of  the  RF 
currents  and  JlQ  and  combining  the  linearized  Maxwell's  equations  (10)  - 
(11)  allows  one  to  write  a  single  differential  equation  for  $=rElQ: 


+  -  0 

where 


(14) 
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(15) 


A  -  i  (1  +  C/(%)) 


e  ■  St<S  -  4>r  !?tv~,D'l«lth»  +  2eVV't(w») 


2  -U-l. 


*o  ° 


2  2 


+a_1d-1(c/Yo)  (^-f  -i2-^2)  +  1 
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(16) 


o>2r2  .2  j.  n-lr/«2  2,  2. 
— 5  -l  +  D  C(fi  r  /c  ) 


(17) 


d  -  pq  -  n2/J  -  i -s/y2  +  y2h2-  n2/u2 


(18) 


and  other  symbols  have  been  defined  previously.  [Equation  (14)  Is  Identical 
to  Eq.  (4)  of  Ref.  28,  except  for  a  few  differences  In  notation.] 

The  other  RF  fields  are  given  In  terms  of  $  by 

Elr  "  lrA  t(F  +  *D  «V0/C  >♦'  2  C»-1(l-Ht)*J  (19) 

c 

*l2  “  -  ^  A_1[(«  +  CD"1!!)*'  -  A  UoCD_1(l4h)t]  (20) 

c 

and  the  perturbed  velocities  may  be  easily  recovered  using  Eqs.  (19)  and  (20) 
In  Eqs.  (13a, b);  the  perturbed  density  Is  obtained  from  ^“-(e^/e)^^. 

Equation  (14)  must  be  solved  subject  to  $(a)  -  $(b)  -  0  (If  the  walls  are 
perfect  conductors)  In  order  to  obtain  the  eigenvalue,  u>.  Hereafter  we  shall 
assume  that  the  electron  layer  is  restricted  to  the  annular  region  rt<r<r2 
[Fig.  1],  In  this  case  it  Is  convenient  to  formulate  the  problem  using  the 
explicit  vacuum  solutions  at  r»r1  and  r-r2  and  to  match  the  appropriate 
logarithmic  derivatives,  as  is  often  done  In  microwave  tube  theory,  in  order 
to  isolate  the  effect  of  wall  boundary  conditions.  Defining,  then,  the 


normalized  admittances  evaluated  just  outside  of  the  layer,  we  have 
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We  stress  that  the  normalized  wave  admittances  b+  and  b_  depend  only  on  the 
waveguide  geometry  exterior  to  the  electron  beam.  They  are  independent  of  the 
beam  or  Its  dynamics.  b+  and  b_  are  evaluated  explicitly  in  Appendix  B  for 
some  practical  cases  of  interest. 

Let  us  examine  Eq.  (1A)  more  closely.  In  the  complex  r  plane  Eq.  (1A) 
has  singularities  at  points  where 


0)  -  Jtdl  *  0 
o 

(22) 

d  +  lh20  -  o 

(23) 

A  -  0. 

(2A) 

The  first  of  these  clearly  represents  a  match  between  the  mode  frequency  and  a 
harmonic  of  the  particle  "cyclotron"  frequency.  Such  a  match  is  present  in 
the  negative  mass  instability,  cyclotron  maser  Instability,  and  diocotron 
Instability  all  of  which  are  described  by  Eq.  (14).  Indeed,  these  may  all  be 
considered  to  be  the  same  Instability  In  this  sense,  though  the  Individual 
names  are  still  useful.  The  negative  mass  instability  is  fundamentally  a 
rotational  effect.  The  classical  explanation  attributes  It  to  the  decrease  of 
circulation  frequency  with  particle  energy  leading  to  growth  of  azimuthal 
bunches.  The  negative  mass  Instability  operates  without  regard  for  the  cavity 


nodes  of  the  vacuum  chamber;  the  fields  are  not  well  approximated  by  the 
cavity  mode.  The  Instability  is  strongly  enhanced  however  if  Re(u>)  happens  to 
be  a  cavity  mode,  as  first  pointed  out  by  Briggs  and  Neil.10  In  recent  years 
this  synchronous  case  has  been  given  its  own  name,  the  cyclotron  maser 
instability,  mainly  in  the  literature  of  gyrotron  research.17*1® 

The  diocotron  instability,  originally  studied  in  connection  with  the 
development  of  the  first  cross  field  microwave  devices,  is  a  "residual" 
instability  in  this  context.  It  is  what  survives  in  the  non-relativistic 
and/or  planar  limits.  The  perturbed  electric  fields  are  basically 
electrostatic  in  nature  and  are  strongly  localized  at  the  position  of  the 
layer.  The  equation  governing  the  diocotron  instability1***®*7*29  may  be 
recovered  formally  from  Eq.  (14)  by  taking  c-*-».  It  Is  discussed  further  in 
Section  IV,  below. 

A  mode  satisfying  Eq.  (22)  is  sometimes  called  the  longitudinal  mode 
since  the  major  effect  on  the  beam  is  azimuthal  bunching.  Equation  (23),  on 
the  other  hand,  describes  what  has  been  called  the  transverse  mode  which 
involves  little  bunching,  but  significant  transverse  (radial)  motion  of  the 
beam.  If  h  is  small,  Eq.  (23)  becomes 

“  "  "  *  V  (25> 

This  transverse  mode  has  been  invoked  to  explain  the  operation  of  the 
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peniotron  J  when  u>  corresponds  to  a  cavity  mode.  In  the  nonsynchronous  case 
the  transverse  mode  is  stable  In  the  absence  of  resistive  walls,  as  pointed 
out  in  ref.  10. 

Finally  we  note  that  Eq.  (24)  may  be  loosely  associated  with  an 
electromagnetic  mode.  In  vacuum  only  (24)  survives  as  a  singularity  but  in 


the  vacuum 


fact  It  may  be  shown  to  be  only  an  "apparent”  singularity,  0 
eigenfunctions  are  analytic  in  r,  except  possibly  at  r-0. 

Ve  proceed  to  analyze  Eq.  (14)  to  uncover  the  parametric  dependences  of 
the  growth  rates  of  the  unstable  modes*  In  order  to  make  progress 
analytically  we  consider  the  case  of  a  beam  of  uniform  density  and  of 
thickness  x  =  r^  -  r^  which  is  much  less  than  its  average  radius 
R  =  (r2  +  II  the  beam  is  sufficiently  thin  it  is  possible  to  consider 

a  Taylor  series  solution  about  r-R  of  Eq.  (14),  taking  care  to  check  at  the 
end  that  the  singularities  are  sufficiently  far  away  in  the  complex  r  plane, 
so  as  not  to  disturb  the  series  convergence;^  carrying  out  this  program  we 
find  the  dispersion  relation 


-  <=.  +  i  <V  G.K^i  +  I)  +  c+cj 

R  A  o 

o 

C 

-  M)2<g+  -  ■ 0 

o 


correct  to  order  (t/R)  .  In  Eq.  (26)  G_,  G+  are  the  radial  logarithmic 
derivatives  of  4  evaluated  just  inside  the  layer  at  r-r^,  r*^  respectively, 
obtained  by  integrating  Eq.  (14)  across  the  beam  edges: 


:*'r 

r2,lMr2,l; 


where 


q(r)  -  «  (£  -  4)  ii l±±l 
1  n  V(I)  y%'  AD 
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X(r)  *  Aq  +  (r-R)A1  +  •••  (29) 

2(r)  *  c0  +  *  V 

near  r-R. 

Expanding  G±  to  order  t/R,  Eq.  (26)  then  reduces  to 


j<b+  +  b_)  +  i  [-Rq'  -  j-(|b+-  q)(|b_+  q)  +  Cj-0 


(30) 


where  now  q  and  q7  are  understood  to  be  evaluated  at  r=R  and  where  in  Eq.  (30) 
we  have  (temporarily)  kept  terms  only  to  0(t/R).  No  assumption  has  been  made 
about  the  E-layer  density  or  current;  we  have  assumed  only  that  the  beam  is 
"thin". 

Let  us  first  consider  the  longitudinal  mode,  that  is,  we  look  for  a  root 
2  2 

of  Eq.  (30)  with  |fl  |«<i>o.  the  terms  within  the  brackets  in  Eq.  (30)  only 
t  2  ~2 

the  q  and  q*  terms  behave  as  fl  for  small  R.  Keeping  only  such  terms,  some 
algebra  gives  Ri=Im(n): 


a.  = 


.  T  2 

1  r°P 


(*o  + 


2h) 


?V+  b->  <1+T2h2) 


+  0Ut/R)‘ 


(31) 


which,  for  h  «  0,  is  the  classical  negative  mass  dispersion  relation  for  the 
Astron  configuration.10  (The  0(S,t/R)2  term  is  displayed  fully  in  Eq.  (36), 
below.)  We  remark  that  Eq.  (31),  including  the  second  order  terms,  can  be 
obtained  from  the  much  simplified  form  of  Eq.  (26), 


G+-  G_  +  xG+G_  -  0 


(26') 


where  terms  behaving  as  G  have  been  dropped  from  Eq.  (26). 
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Ignoring  the  0(lt/R)  term  in  Eq.  (31)  for  a  moment,  we  may  Interpret  the 

dispersion  relationship  (31)  as  follows.  As  we  have  seen,  the  factor 

(b+  +  b_)  represents  the  effects  of  the  container  structure,  and  the 
2 

factor  (T/R)b>p  Is  proportional  to  the  beam  current.  These  two  terms  always 
appear  in  any  dispersion  relation  involving  an  electron  beam.  All  dynawlcal 
effects,  including  those  due  to  self  fields,  are  contained  In  the  factor 


M  =  - 


+  2h) 

2  2 

(i  +  yV) 


for  a  sufficiently  thin  beam.  Some  insight  into  the  meaning  of  this  factor 

may  be  had  If  we  note  that  for  small  equilibrium  densities,  £,  It  is 

2 

proportional  to  dwo/de  where  e  *  me  is  the  total  energy  of  a  particle 

in  the  equilibrium: 

dui  -  w  +  2h  +  £/y* 

O  Or  O  Oi  .... 


Or  o _ 

W  +  tV" 


At  least  when  £  is  small  then  we  find  agreement  with  the  classical  explanation 

of  the  negative  mass  effect:  A  rotating  beam  is  unstable  (stable)  if  its 

equilibrium  rotation  frequency  is  a  decreasing  (Increasing)  function  of  its 

(total)  energy,  that  is,  M  is  proportional  to  the  effective  azimuthal  inertia 

(mass).  We  stress  that  for  finite  5  this  interpretation  begins  to  break  down; 

Eq.  (31)  was  derived  assuming  £  was  finite  yet  £  does  not  appear  in  the  factor 

M  but  clearly  does  appear  in  Eq.  (33). 

M,  though  simple  in  form,  has  many  interesting  features.  Let  us  consider 

its  dependence  on  h,  illustrated  in  Fig.  3.  Perhaps  the  most  interesting 

2 

property  of  M  is  that  it  experiences  a  change  in  sign  at  h  *  -  B  /2.  The 

o 

possibility  for  using  this  property  of  equilibria  supported  by  a  radially 
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Inward  electric  field  in  addition  to  a  magnetic  field  in  order  to  suppress  the 

negative  mass  instability  has  been  discussed  elsewhere. We  note  that  since 

the  factor  M  is  independent  of  the  beam  current  or  of  the  waveguide  dimensions 

and  represents  a  purely  dynamical  quantity,  this  stabilisation  of  the  negative 

mass  instability  by  a  radial  electric  field  is  expected  to  be  valid  even  in 

toroidal  geometry,  even  for  high  current,  very  cold  beams.  (The  method  is 

limited  to  use  for  moderately  low  energy  beams,  in  the  MeV  range  in  practice, 

since  the  applied  electric  field  required  to  change  the  sign  of  M  from 

3 

negative  to  positive  is  proportional  to  y  »  as  follows  from  the  definition  of 

h,  Eq.  (5).)  Note  that  there  is  a  "most  stable"  configuration  at  h— 1,  at 

which  point  certain  singular  parts  of  C  in  Eq.  (14)  vanish  identically. 

Perhaps  equally  interesting  for  microwave  generation  applications,  is  the 

occurrence  of  a  "most  unstable"  point  in  Fig.  3  at  h  =  l/Y2.  If  we  recall 

o 

2 

that  yQh  is  the  ratio  of  the  radial  electric  force  to  the  centrifugal  force  we 

2 

observe  that  the  choice  yQh  *  1  describes  a  configuration  in  which  the 

equilibrium  is  supported  solely  by  an  electric  field  (a  Type  III  equilibrium, 

in  the  language  of  Section  II).  This  result  suggests  that  for  a  given  beam 

energy,  a  microwave  source  such  as  the  orbitron^  in  which  an  annular  beam 

circles  a  positively  charged  wire  with  no  applied  magnetic  field,  might  have 

some  advantage  over  more  conventional  devices  like  magnetrons  (h<0),  inverted 

magnetrons  (h»l/yQ),  or  gyrotrons(h«*0) .  This  finding  is  novel.  We  remark 

o 

that  this  negative  peak  in  M  is  sharpest  at  low  energies  (small  6  )• 

o 

2  2 

In  the  non-relativistic  limit,  B  -*•  0,  y  -►  1,  the  sign  of  M  is  determined 

o  o 

directly  by  the  sign  of  h,  which  remains  finite  in  this  limit.  Usually  we 
are  accustomed  to  thinking  that  the  negative  mass  instability  should  vanish  in 
the  non-relativistic  limit  when  there  is  no  gradient  in  the  magnetic  field. 


Here  we  see  however  that  an  equilibrium  electric  field  can  affect  the  sign  of 
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the  azimuthal  inertia  just  as  a  gradient  in  B  can.  If  we  define  an  equivalent 
field  index,  neq,  by 


then 


eq 


2  <Y?-Uh 

h[—  -P— -1 

(1  +  h) 


(34) 


(35) 


which  vanishes  when  h  does. 

One  final  comment  about  h  must  be  made.  The  electric  field  appearing  in 
the  definition  of  h,  Eq.  (5),  is  due  both  to  the  electric  charge  of  the  layer 
and  to  any  externally  applied  bias  potential.  When  the  contribution  of  the 
electric  charge  dominates  there  is  some  question  as  to  what  value  of  h  to  use 
In  Eq.  (31)  when  calculating  growth  rates  since  Eq  and  therefore  h  will  change 
6ign  somewhere  within  the  layer.  For  this  reason,  in  all  numerical  examples 
considered  below  we  will  assume  that  h,  when  it  Is  non-negligible,  is 
dominated  by  the  contribution  of  an  externally  maintained  bias  potential. 

The  analytical  results  presented  up  to  this  point  are  subject  to  test  by 

numerical  integration  of  Eq.  (14)  subject  to  suitable  boundary  conditions.  We 

have  written  a  program  to  carry  out  this  task.  Given  a  density  profile,  the 

electrostatic  potential  difference  and  the  total  magnetic  flux  contained 

between  rsa  and  b,  the  program  calculates  all  equilibrium  quantities  then, 

given  the  mode  number  £,  locates  an  eigenvalue  u.  In  all  examples  below  the 

density  profile  is  parabolic;  we  expect  good  agreement  with  Eq.  (31),  derived 

2  2  2 

using  a  flat  profile,  if  we  identify  (t/R)  with  2(v/yq)(c  /R  )  where  v  is 
Budker's  parameter,  the  number  of  particles  per  unit  length  times  the 
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classical  electron  radius.  All  examples  use  a  perfectly  conducting  wall; 
effects  of  a  small  wall  resistivity,  which  may  be  very  important  under  some 
conditions,  are  discussed  briefly  toward  the  end  of  Appendix  B. 

Though  we  shall  eventually  test  Eq.  (31)  for  all  parameter  dependences, 
let  us  first  check  the  interesting  dependence  on  h,  discussed  above,  against  a 
numerical  solution  of  Eq.  (14).  Figure  4  presents  some  typical  results.  Here 
we  have  illustrated,  for  two  different  values  of  v/yq  a  comparison  of 
the  1=1  mode  growth  rates  for  a  thin  (t/R  =.02)  beam  as  predicted  by  Eq. 

(31)  (dashed  line)  and  by  a  solution  of  Eq.  (14)  (solid  line).  (To  be 
precise,  what  is  plotted  on  the  numerical  solution  curve  is  Im(ci})/u>  (r^).) 

The  transition  from  stable  to  unstable  behavior  occurs  at  h  =  -  B0/2 

(-  -  0.278  for  this  case)  independent  of  £.  This  is  an  important  feature  of 

OQ 

Eqs.  (31)  and  (32)  which  has  been  used  to  argue  0  that  the  stabilization 
2 

condition,  h<-BQ/2,  for  the  negative  mass  instability  is  independent  of  the 

beam  current.  Note  that  this  statement  cannot  be  made  on  the  basis  of  single 

particle  orbit  theory  alone  [cf.  Eq.  (33)];  the  stabilization  condition  is 

obtained  from  collective  mode  considerations  including  self  field  effects. 

Actually  a  small  amount  of  residual  growth  remains  for  h  less  than  but  very 
2 

close  to  -  8  / 2,  a  feature  discussed  in  more  detail  in  the  following 
o 

section.  Agreement  between  Eq.  (31)  and  the  numerical  solution  is  best  at 

small  growth  rates  which  is  reasonable  if  we  recall  that  terns  behaving 

-1  -2 
as  fi  were  neglected  in  favor  of  terms  behaving  as  when  Eq.  (31)  was 

derived.  Equation  (31)  consistently  gives  slightly  too  high  a  growth  rate 

(i.e.  it  is  pessimistic),  a  feature  which  will  be  shown  to  persist  when 

variations  of  other  parameters  are  considered.  Notice  that  the  numerical 

solution  confirms  the  existence  of  a  peak  in  the  growth  rate  for  h  =  1/yo»  as 

predicted  by  Eq.  (31),  again  independent  of  £. 


The  significant  parameters  upon  which  the  small  signal  behavior  of  the 
layer  depends  include,  beside  the  externally  applied  bias  fields,  the  geometry 
factors  a,  b,  r^,  r2»  the  current  (v),  the  beam  energy  yq,  and  the  mode 
number  £.  We  proceed  to  consider  the  effect  of  each  of  these  separately  on 
the  growth  rate  of  the  longitudinal  mode,  u  -  to>o,  for  the  specific  case  of  an 
Astron-like  configuration,  h  «  0. 

Unless  otherwise  stated  all  parameters  in  the  cases  considered  below  will 

take  the  following  nominal  "base  case"  values:  a  -  0.5m,  b  -  2.2m,  rj  - 

0.99m,  r2  ■  1.01m,  total  electrostatic  potential  difference  between  inner  and 

outer  walls  -  0.,  total  flux  between  inner  and  outer  walls  -  *(b2  -  a2}  B 

v  ’  oo 

-4 

where  BOQ  =  48.2x10  T  is  the  field  required  to  hold  a  single  particle  at 

R  -  1m  with  y  “  3.  The  radial  density  profile  is  always  taken  to  be  parabolic 

and  symmetric  about  r«R=lm  with  specified  peak  value;  the  base  case  value  is 
7  _3 

5x10  /cc  which  gives  v/yq  *  3.94x10  .  The  base  case  azimuthal  mode  number 
£  Is  1,  for  which  b+  +  b_  -  2.50  for  u>  -  u>o.  From  Eq.  (31)  the  normalized 
growth  rate  for  the  base  case  is  5.6Z;  the  numerical  solution  of  the 
eigenvalue  problem  gives  5.35Z. 

Figure  5  illustrates  the  comparison  of  growth  rates  as  calculated  by  Eq. 

(31)  and  by  a  numerical  solution  of  the  eigenvalue  problem  for  a  range  of 

1/2 

currents.  Over  the  wide  range  considered,  the  v  scaling  predicted  by  Eq. 
(31)  is  shown  to  hold  up  extremely  well  up  to  values  of  v/yq  of  a  few 
percent.  Similar  excellent  agreement  is  generally  found  for  variations  in 
layer  thickness  (Fig.  6),  outer  wall  position  (Fig.  7),  inner  wall  position 
(Fig.  8),  particle  energy  (Fig.  9),  and  azimuthal  mode  number  (Fig.  10).  Some 
remarks  on  each  case  follow. 

In  the  case  of  varying  layer  thickness  with  the  maximum  density  fixed 
(Fig.  6)  two  effects  are  competing;  these  are  the  basic  v1^2  dependence  of  the 


growth  rate  versus  the  stabilizing  effect  of  finite  thickness30  (effectively, 

finite  frequency  spread).  The  finite  thickness  effect  is  clearly  second  order 

in  (t/R)  and  is  not  shown  explicitly  in  Eq.  (31)  but  it  included  in  all 

plotted  data  of  Figs.  5  -  10;  the  second  order  term  is  given  explicitly  in  the 

following  section.  In  the  examples  we  have  studied,  finite  thickness  effects 

have  been  small  and  have  not  been  effective  in  stabilizing  the  instability. 

In  Fig.  7  the  effect  of  outer  wall  position  is  Illustrated.  As  the  wall 

is  moved  in  from  its  base  case  location  at  2.2m  the  growth  rate  is  observed  to 

increase  dramatically  for  a  while,  then  to  fall  off.  The  reason  for  the 

increase  is  the  approach  to  synchronism  of  particle  motion  with  a  cavity 

vacuum  mode,  that  Is,  u  «  Lwq  «  a>v  where  mv  is  a  solution  to  b+  +  b_  -  0; 

under  this  synchronous  condition  the  cavity  mode  fields  act  to  enhance  those 

established  by  the  dynamical  charge  bunching  due  to  the  negative  mass  effect. 

This  synchronous  case  has  been  given  its  own  name,  the  cyclotron  maser 

instability,  and  Is  put  to  enormous  practical  use  in  the  gyrotron  family  of 

microwave  devices.17  While  in  Figs.  5  -  10  we  have  consistently  evaluated 

b+  and  b.  at  u  •  Iu>o,  near  a  zero  of  b+  +  b_  this  Is  clearly  inadequate  and 

Eq.  (31)  should  be  solved  as  a  cubic  polynomial.  Empirically  we  find  that 

evaluating  b+  +  b.  at  u  *  iuo  Is  adequate  when  b+  +  b  >  1 . 

As  the  outer  wall  is  moved  further  inward,  past  the  synchronous  point  the 

growth  rate  in  Fig.  7  drops  as  b+  +  b_  changes  sign.  This  drop  is 

attributable  to  a  "shorting  out"  of  the  azimuthal  field  E.  as  the  wall 

19 

approaches  the  edge  of  the  layer.  An  identical  phenomenon  is  seen  as  the 

inner  wall  is  moved  outward.  Fig.  8.  Use  of  an  inductive  impedance  (b+  +  b_ 

<0)  to  stabilize  the  negative  mass  instability  has  been  proposed  by  Briggs  and 
32 

Neil.  We  remark  that  in  Fig.  8  no  synchronous  case  is  encountered  for  the 
parameters  considered,  as  the  inner  wall  is  moved. 
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The  decrease  of  growth  rate  with  increasing  kinetic  energy  ie  documented 
in  Fig.  9.  The  basic  reason  for  this  decrease  is  just  the  relativistic  mass 
Increase:  the  azimuthal  response  of  a  particle  to  the  perturbed  field  E  is 

In 

reduced  by  a  factor  ,  for  large  y0*  No  synchronous  cases  are  encountered 
over  the  range  of  y  considered. 

In  Fig.  10  we  have  plotted  growth  rate  versus  azimuthal  mode  number  l. 
Agreement  is  good  between  the  dispersion  relation,  Eq.  (31),  and  the  numerical 
solution;  near  synchronous  cases  occur  for  g,  -  5  and  9.  Though  predicted 
growth  rates  are  rather  large,  for  the  high  i  modes,  these  should  in  practice 
be  subject  to  stabilization  by  the  effects  of  finite  betatron  oscillations 
which  we  have  neglected  in  this  treatment.  In  any  event  we  expect  the 
dispersion  relation,  Eq.  (31),  to  begin  to  break  down  for  short  wavelengths, 
i.e.  £t/R  ~  0(1),  or  azimuthal  wavelengths  on  the  order  of  the  layer 
thickness.  A  WKB  treatment  of  Eq.  (14)  carried  out  for  this  case  yields  an 
eigenvalue  condition,  the  numerical  solution  to  which  would  appear  to  require 
more  effort  than  a  direct  numerical  solution  of  Eq.  (14)  Itself. 

Finally,  we  report  that  in  the  initial  phase  of  our  investigation,  we 
attempted  to  derive  the  dispersion  relationship  by  constructing  a  quadratic 
form  for  the  differential  equation  (14),  under  the  assumption  that  $  -  rE.  is 
approximately  constant  across  the  E-layer.  The  latter  assumption  has  been 
widely  used  and  has  been  considered  valid^”**»^»^  for  a  thin  relativistic 
electron  beam.  However,  this  line  of  investigation  led  to  an  incorrect 
dispersion  relationship  If  up/mo  ~  0(1),  even  to  the  lowest  order  in  t/R.  In 
other  words,  to  account  correctly  for  the  DC  self  fields  in  the  present 
Eulerian  description,  the  tangential  AC  electric  field  should  not  be  assumed 
constant  across  the  E-layer,  regardless  of  the  thickness.  On  the  other  hand, 
our  dispersion  relationship  Eq.  (31)  correctly  accounts  for  the  self  fields, 


and  is  valid  for  arbitrary  beam  energy,  and  arbitrary  combination  of  EQ  and 
Bq>  as  already  stressed. 

IV.  Residual  Growth:  The  Diocotron  Instability 

Sufficiently  close  to  the  "zero  mass"  points  h  -  -  g  / 2,  h+  ± 

o 

Eq.  (31)  begins  to  be  dominated  by  the  neglected  terms  of  order  (tVR)2.  We 
have  already  observed  this  phenomenon  near  the  point  h  *  -  6/2  in  Fig  4.  In 
this  section  we  discuss  the  point  h  -*■  ±  «*,  which  extreme  is  reached  in  the 
planar  limit:  r  -*■  •,  t  •  ,  Z/r  =  ky  finite.  We  remark  that  the  vanishing  of 

the  negative  mass  growth  rate  in  the  planar  limit,  a  feature  of  Eq.  (31)  and 
clearly  expected  on  physical  grounds,  has  not  previously  been  demonstrated 
analytically,  to  the  authors'  knowledge. 

The  second  order  term  in  Eq.  (31)  is  given,  in  complete  generality, 
by  (At/R)2oj2A  where 


4Y^(l+h)2(l+Y2h~) 


1  T  U2l»o  +  2h  *  <1+h>2J 


+2^o  (Bo  +  2h)2  *  ^Cl^yXB2  +  2h)2  ] 


b  -  b 


_  (1  +  h)  (jf _ _^.i2 

l2  2. .  .  2.  2,  V  +  b  ti 

V1  +  V1  5  + 


(36) 


o 

As  it  stands  A  includes  both  the  diocotron  instability  (£  terms)  and  the 
finite  thickness  stabilization  effect^0  ((°  term)  referred  to  above.  The  last 
term  of  Eq.  (36)  contains  wall  boundary  effects;  specifically  we  may  recover 
the  stabilization  at  the  diocotron  mode  due  to  contact  of  the  layer  with 
perfectly  conducting  walls, B  as  well  as  the  destabilizing  effects  of  finite 
wall  conductivity . ^  ^  If  we  assume  that  the  walls  are  many  wavelengths  away 
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■  1. 


and  that  the  fields  are  electrostatic  then,  in  the  planar  limit,  b^ 

3  2 

Taking,  then,  the  planar  limit,  hw  -  eB  /(my  )  =  -  u  /y  ,  only  the  first 

w  u  O  CO 

2  L 

£  term  in  Eq.  (36)  survives  and  we  recover  the  classical  result*  for  the 
growth  rate  of  the  dlocotron  mode 


Im(w)  -  |k  tv*/2| 

y  o 


(37) 


p  2  2 

where  the  velocity  shear  v  »  u  /(w  y  ).  Of  course  Eq.  (37)  may  be  obtained 

O  p  CO  *  * 


by  much  simpler  methods  than  those  employed  here;  our  point  is  only  that  it  is 


recoverable  from  the  present  formalism.  Note  that  the  dependence  on  the  line 


density,  v,  of  the  dlocotron  growth  rate  is  just  v  ;  for  the  (non-synchronous) 

1/2. 


negative  mass  Instability  the  dependence  is  v  ;  for  the  synchronous  case  it 
1/3 


is  v 


The  relationship  between  the  dlocotron  and  negative  mass  instabilities 

has  been  discussed  by  Neil  and  Heckrotte^®  and  by  Lau.*®  Mostrom  and  Jones^ 

have  recently  examined  the  electrostatic  case,  including  the  effects  of  shear 

29 

in  vz.  Davidson  and  Tsang  *  have  reported  analytical  and  numerical  results  in 
cylindrical  geometry. 


V.  The  Transverse  Mode 

An  electron  moving  in  a  field  satisfying  Eq.  (23)  where  Re(u)  is  a  vacuum 
guide  mode  may  be  shown  to  be  acted  upon  by  a  nearly  constant  electric  field, 
when  the  motion  is  averaged  over  its  gyro-orbit.  The  particle  therefore 
experiences  an  t  x  t  drift  transversely,  toward  the  wall,  which  motion  brings 
the  particle  to  experience  yet  a  stronger  electric  field.  A  net  transfer  of 
energy  from  (to)  the  particles  may  be  shown  to  result  for  the  -  (+)  sign 
resonance  of  Eq.  (25).  This  mechanism  of  wave  growth  has  been  used  to  explain 


2  5 

the  operation  of  the  peniotron  oscillator. 

In  the  planar  Unit  of  the  previous  section,  the  resonance  condition  for 
the  transverse  mode,  Eq.  (23),  becomes  y0(u>-  kyvc)  -  ±  uc,  which  is  the  well- 
studied  mode  of  planar  magnetron  tube  theory .1_7*33,34  factor  [£q. 

(5)1,  thus  again  appears,  as  it  did  In  Eq.  (31),  as  a  measure  of  the 
"planarity*'  of  the  configuration. 

A  dispersion  relation  for  the  transverse  mode  may  also  be  obtained  from 
Eq.  (30).  Keeping  the  most  important  terms  we  have,  approximately 


+  bJ  "  I  7T  <jV  q)(K  +  <*>  “  0  (38> 


where,  near  the  zeroes  of  A 


-  (a2-  oV) 

*o  *  ,2  2,,  2  »~2  <39> 

1  V1 '  eois;) 

«-£■  "VI  <«> 

o 

2  2  2 

where  a  *  1  +  y  h  .  Equation  (38),  using  Eqs.  (39)  and  (40),  agrees  with 
Eq.  (39)  of  Briggs  and  Neil10  for  h  ■  0,  to  leading  order  in  £• 

Clearly  there  are  no  unstable  roots  of  Eq.  (38)  near  the  (simple)  zeroes 
of  AQ  unless  either  one  of  these  nearly  coincides  with  a  zero  of  b+  +  b_ 

(a  guide  mode  cutoff  frequency)  o£  b+  or  b_  contains  an  imaginary  part,  due  to 
finite  wall  resistivity  for  example  (see  Appendix  B).  Thus  the  transverse 
mode  (for  small  h)  depends  crucially  on  the  interaction  of  the  electrons  with 
their  external  surroundings,  unlike  the  negative  mass  and  diocotron 
instabilities,  the  mechanisms  of  which  operate  in  a  manner  that  is  insensitive 
to  boundary  conditions  on  the  fields  at  distant  walls.  In  the  synchronous 


case,  with  perfectly  conducting  walls  the  predicted  scaling  of  the  growth  rate 


of  the  transverse  node  with  v  is  v 


The  dependence  of  the  real  part  of  the 


frequency  of  the  mode  on  an  externally  applied  radial  DC  electric  field  nay 
prove  useful  in  some  circumstances. 

The  penlotron  interaction  is  essentially  non-relativistic.  It  relies 
heavily  on  the  spatial  Inhomogeneity  of  the  perturbed  wave  fields.  In 
contrast,  for  longitudinal  modes,  the  spatial  inhomogeneity  of  the  unperturbed 
motion  (l.e.  shear)  is  far  more  important.  Both  the  transverse  and 
longitudinal  modes  can  be  used  to  convert  the  rotational  energy  of  the 
electrons  to  rf  waves  efficiently,  however. 

In  the  planar  limit  Eq.  (38)  continues  to  predict  stability  in  the 
absence  of  a  resistive  wall  yet  it  is  well-known*-^  that  inclusion  of  a 
resonant  layer  satisfying  Y0(u»  -  kyvQ)  *  ±  a>c  in  the  beam  leads  to  wave 
growth.  The  resolution  of  this  contradiction  lies  in  the  failure  of  the 
Taylor  series  solution  to  Eq.  (14),  from  which  Eq.  (38)  was  obtained;  the 
resulting  growth  rate  is  non-analytic  in  £  and  one  must  resort  to  numerical 
or  other  methods^  to  solve  the  eigenvalue  problem. 
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VI.  Summary 

In  this  paper  ve  have  attempted  a  general  treatment  of  the  linear 
stability  problem  for  laminar  electron  flow  In  cylindrical  geometry.  The 
basic  equilibrium  state  has  been  taken  to  be  maintained  by  radial  electric 
and/or  axial  magnetic  fields  (Eq.  (1)).  No  azimuthal  magnetic  field  nor  any 
axial  electron  motion  have  been  included  In  the  equilibrium  state.  The  linear 
stability  problem  for  azimuthal  and  radial  perturbations  has  been  formulated 
exactly,  fully  relativistically  and  fully  electromagnetically ,  including  all 
effects  of  self  fields.  The  stability  problem  reduces  to  an  eigenvalue 
problem  for  the  frequency  u,  given  the  azimuthal  mode  number,  £  (Eq.  (14), 
with  associated  boundary  conditions). 

Our  efforts  have  been  focused  on  the  longitudinal  mode  (Eq.  (22)),  for 
thin  beams,  which  is  of  considerable  importance  in  accelerator  and  microwave 
device  research.  We  have  obtained,  and  favorably  compared  to  a  numerical 
solution  of  the  eigenvalue  problem,  a  dispersion  relation  In  the  thin  beam 
limit  (Eq.  (31))  which  applies  in  complete  generality  to  the  longitudinal  mode 
and  which  reproduces  all  classical  results  in  appropriate  limits.  Some 
interesting  differences  among  equilibria  regarding  the  negative  mass 
instability  have  been  pointed  out;  namely  ve  have  found  a  simple  way  either  to. 
stabilize  or  to  maximize  the  growth  of  this  mode.  This  finding  might  have 
practical  consequences  in  accelerator  or  microwave  tube  design. 

The  longitudinal  mode,  a  -  luo»0,  encompasses  the  negative  mass,  electron 
cyclotron  maser,  and  diocotron  instabilities.  The  negative  mass  and  electron 
cyclotron  maser  effects  are  unique  to  cylindrical  geometry;  they  are 
fundamentally  relativistic  in  nature  when  the  motion  Is  supported  solely  by  a 
magnetic  field.  They  are  even  more  pronounced,  especially  for  low  energy 
beams,  if  the  equilibrium  rotation  Is  supported  solely  by  a  radial  electric 
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field.  In  planar  geometry  both  of  these  instabilities  are  absent  and  only  the 
residual  dlocotron  instability  remains  which  itself  may  be  stabilized  by 
placing  the  layer  in  contact  with  a  conducting  wall.^ 

The  transverse  mode,  w-£o>0  ■  i  otao  (a  is  defined  following  Eq.  (40))  has 
been  used  to  explain  the  operation  of  the  peniotron  device.  When  the  geometry 
is  cylindrical  (|yoh|«l)  this  mode  is  stable  unless  the  electron  motion  is 
synchronous  with  a  cavity  mode  and/or  resistive  walls  are  present.  In  the 
planar  limit  o^0*^c/'T0  ®nd  we  identify  this  mode  as  the  Doppler  shifted 
cyclotron  resonance  considered  by  Buneman^'5  and  others1, 6,^»33  in  stU(jieB  ef 
magnetron  operation.  This  mode  is  the  dominant  unstable  mode  for  planar,  high 
density  laminar  flow. 

Finally  we  remark  that  the  singularities  defined  by  Eq.  (24),  which  we 

have  not  examined  here,  may  be  worth  6ome  additional  study;  however  we  note 

30 

that  in  both  the  vacuum  case,  £  ■  0,  and  in  the  case  of  planar  Brillouin 
flow, 1,2,5,33  tu*  -  us\  the  singularity,  Eq.  (24),  is  removable. 
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Appendix  A 


Technique  for  Numerical  Solution  of 
Equilibrium  Problem 

In  this  Appendix  we  describe  a  simple  technique  based  on  Newton's  method 
which  we  have  used  to  solve  a  certain  eigenvalue  problem  associated  with  the 
calculation  of  laminar  E-layer  equilibria.  Given  the  electric  field  which  is 
trivially  solved  for,  having  specified  the  density  profile  and  potential  on 
one  wall,  we  must  find  the  momentum  profile,  u(r)  and  the  field  BQ(r)  subject 
to  the  constraint  of  specified  total  flux  Eq.  (7).  The  total  flux  might  be 
specified  in  an  experiment  in  which  a  beam  is  injected  and  contained  in  a 
chamber  for  less  than  a  magnetic  diffusion  time. 

Let  us  cast  the  problem  in  terms  of  the  dimensionless  fields,  a£  and  otg 
defined  in  the  text  following  Eq.  (6).  The  problem  then  is  to  find  u(r)  and 
dg(r)  subject  to 

b 

/  dr  Og  =  F  =  specified  constant  (A-l) 

a 


where  u  and  ag  satisfy 


1-5/Y2  "  Y2h 
ur  o  ’o  t 

u  1  +  h  J 

4  V'  ■  \ u 

c 


(A-2) 


(A-3) 


u2  -  <xe(u2+1)1/2+  aRu 


(A-4) 


for  r^r^rj*  Equation  (A-2)  guarantees  that  if  u(rj)  satisfies  Eq.  (A-4)  then 
u(r)  will  do  so  for  all  r.  The  algorithm  proceeds  as  follows:  An  initial 
guess  is  made  for  a_(r.)  using  the  value  of  the  externally  applied  B  field, 

D  1 
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say.  Using  the  known  value  of  (^(r^)  the  roots  (two,  in  general)  of  Eq.  (A-4) 
are  found  and  the  one  corresponding  to  the  equilibrium  of  Interest  Is  chosen. 
Equations  (A-2)  and  (A-3)  are  then  Integrated  and  the  difference 


/  dr  aB-  F  =  D* 

a 


(A-5) 


Is  calculated,  ag(r^)  Is  then  adjusted  according  to 


°B,n+l  "  aB,n  b  3a-  (r) 

fdr- -  _ 

I  3aB,n<'l> 


(A-6) 


where  the  subscript  n  denotes  an  iteration  number.  The  loop  is  stopped  once 
D*  is  less  than  some  specified  tolerance. 


It  remains  only  to  describe  the  evaluation  of  the  denominator  in  Eq. 

(A-6):  The  dependence  of  aB(r)  on  its  initial  condition  at  ^  Is  found  during 

the  integration  of  Eqs.  (A-2)  and  (A-3)  by  simultaneously  integrating  the 

equations  for  3oB(r)/3ctB(r1)  and  3u(r)/3<*B(rl)  which  are  simply  obtained  by 

explicitly  differentiating  Eqs.  (A-2)  and  (A-3)  with  respect  to  a„(r.).  The 

B  1 

initial  condition  3u(r. )/3aB(r. )  is  obtained  from  Eq.  (A-4). 


Appendix  B 


Evaluation  of  the  TE  Wave  Admittances 
In  this  Appendix  the  normalized  admittances  referred  to  in  the  text  in 
Eqs.  (21a, b)  are  given  for  the  geometry  of  Fig.  1;  the  toroidal  and  planar 
cases  are  also  discussed. 

In  the  vacuum  regions  a<r<rj  and  r2<r<b  the  wave  equation  Eq.  (14)  is 


w  r  .2 

2 

c 


(B-1) 


the  general  solution  to  which  is 


♦  -  x  [C^Cx)  +  C2Y'(x)] 


(B— 2  ) 


where  x  =  ur/c,  and  Y^  are  the  standard  Bessel  functions,  Cj  and  C2  are 

constants,  and  in  this  Appendix  a  prime  will  denote  Note  that  the 

dx 

2  2 

"singularities"  at  x  -  l  in  Eq.  (B-1)  are  only  apparent,  not  real,  i.e.  4>  is 
analytic  at  these  points.  The  other  vacuum  fields  are 


U  d  <♦> 

2  2  dr’ 

x  -l 


(B-3) 


—  A  V 

cl  lr* 


(B-4) 


Using  the  definitions  in  the  text  Eqs.  (21a, b),  it  follows  that  b±  may  be 
generally  written 


C.V(x,)  +  C,+Y  (x0) 


b _ i  ri 2  1 

+  x  l  +  +  J 

2  c.  J'(x,)  +  c,  y:(x,) 


(B-5a) 


(B-5b) 


b  _  »  tCl  W  *  C2  W, 

*1  Cb~J'(Xj)  +  Cj-YJCx,) 


where  “rj  2^c* 


+  + 


The  ratios  /C ^  and  /C^  are  determined  by  the  boundary  conditions 
at  r-b  and  r*a  respectively.  Some  special  cases  of  interest  are: 


1.  Perfectly  conducting  wall  at  r*rw:  (rw=  a  or  b) 


VC2  ’  (xw)/Jt  <V 


(  B— 6> 


2.  Wall  with  (complex)  dielectric  e(ai): 


VC2  "  - 

- 


W  -  mz(y) 


(t^o)Ji(xw>d^lnZ(y) 


]  (B-7) 


where  x.  is  the  surface  impedance,  (p/e)1  ,  co  “  376.70, 
1/2 

y=  xr^,  K=u>  (ep)  ,  Z(y)  =  J^(y)  for  inner  wall, 
H„^^(y)  for  outer  wall  where  H  is  the  Hankel 
function,  c  and  p  refer  to  the  wall  material. 


3.  No  inner  wall: 


C2  "° 


(B-8) 


4.  No  outer  wall: 


Cl+/C2+  -  -  1 


(B-9) 
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b+  =  -  (ky/a)  cot  laCr^a))  (B-ll) 

b_  -  -  (ky/a)  cot  [a(b-r2)J  (B-12) 

where  a  =  (w2/c2-  k2)1/2.  The  planar  limit  of  (B-10)  is 
just  the  electrostatic  limit  of  (B-ll)  and  (B-12),  as 
expected,  in  which  case 


b+  -  coth  [ky(b-r2)J  (B-13) 

b_  -  coth  [ky(r^-a)].  (B-14) 

When  a  wall  is  resistive  the  resulting  dissipation  is  represented  by  an 
imaginary  part  in  the  corresponding  admittance  b±.  (In  Eq.  (B-7)  e(w)-ia/<*> 
for  a  good  conductor  of  conductivity  o>)  That  such  dissipation  can  lead  to 
disruptive  beam  instabilities,  even  for  a  "positive  mass"  beam,  has  been  known 
since  the  pioneering  work  of  Neil  and  Sessler  and  Laslett,  Neil,  and 
Sessler.^  It  is  in  fact  these  resistive  modes,  rather  than  the 
(fundamentally  dynamical)  negative  mass  instability  which  are  thought 
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ultimately  to  limit  the  beam  current  In  cyclic  electron  accelerators. 

If  a  wall  la  not  smooth  but  contains  some  structures  (cavities,  fins, 

etc.)  a  common  practice  Is  to  calculate  an  approximate  value  for  the 

admittance  and  to  use  it  as  a  boundary  condition  for  some  approximation  to  Eq. 

(14)  in  the  text.  In  this  case  however  the  problem  is  not  being  treated  fully 

self  consistently  since  the  equilibrium  of  Section  II  would  not  be  strictly 

correct,  l.e.,  the  correct  equilibrium  would  no  longer  have  azimuthal  symmetry 

(and  would  be  much  harder  to  calculate). 

For  the  case  of  a  perfectly  conducting  wall  we  have  plotted  in  Figs. 

B-1,2,3,4  the  quantity  b+  +  b_  evaluated  for  Xj  -  x2  -  for  various  values 
2 

of  E±  =  me  (yo“1),  assuming  a  thin  layer.  It  is  this  quantity,  b++  b_,  which 
enters  the  dispersion  relation  for  the  longitudinal  mode  discussed  in  the 
text. 

For  a  toroidal  (accelerator)  geometry,  the  dispersion  relation  Eq.  (31) 

Is  expected  to  be  replaced  by 

flj  “  -  e08*2  (2''/Y0)(c2/R2)(mc2/uo)0[,)o/3E)ext  (B-15) 

where  (3<i»o/3e)ext  denotes  the  derivative  of  the  circulation  frequency  with 
respect  to  total  particle  energy  evaluated  as  If  the  particle  were  acted  upon 
only  by  the  external  electric  and  magnetic  fields.  In  Eq.  (B-15),  we  have 
Identified  the  geometric  factor  g  of  refs.  [8,11)  with  l/(le0(b+  +  b  )) ,  as 
pointed  out  in  ref.  10.  Note  that  the  dispersion  relationship  (B-15)  includes 
self  field  effects  and  that  g  is  always  positive  for  toroidal  geometry  (with 
smooth,  perfectly  conducting  walls). 
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Figure  Captions 


Fig.  1.  Model  of  an  E-layer.  The  layer,  of  infinite  extent  in  z  (in  and 
out  of  the  page),  occupies  the  region  rj  <  r  <  ^  between  the  walla 
of  a  coaxial  guide  at  r«a  and  r*b.  The  electrons,  supported  by  an 

A  A 

electric  field  E  (r)r  and  a  magnetic  field  B  (r)z,  move  in  concentric 
o  o 

circles  in  the  equilibrium  state,  either  clockwise  or 
counterclockwise  depending  on  the  equilibrium  type. 

Fig.  2.  The  cig-  aB  plane.  The  number  and  type  of  solutions  to  the 

equation  of  radial  force  balance  (6)  are  shown  in  each  region  and  on 

the  boundaries  between  regions.  A  root  is  labeled  +  or  -  if  it 

corresponds  to  clockwise  or  counterclockwise  rotation  respectively  as 

one  looks  in  the  direction  of  $  . 

o 

Fig.  3.  Plot  of  the  dimensionless  "azimuthal  mass*'  M  versus  the 

dimensionless  equilibrium  electric  field,  h.  The  actual  plot  shown 

is  for  the  case  y  *1.3  but  all  axis  labels  are  expressed  generally 
o 

in  terms  of  y  or  $  . 

o  o 

Fig.  4.  Normalized  growth  rate  for  the  negative  mass  instability  versus  h 
for  the  case  a  *  0.6m,  b  »  2.6m,  R  *  1.0m,  t/R  *  0.02, 

Yo  -  1.5,  l  ■  1.  A  solid  curve  Indicates  data  obtained  from  a 
numerical  solution  of  Eq.  (14);  the  dashed  line  is  a  plot  of  Eq. 

(31).  The  upper  pair  of  curves  is  for  v/yQ  ■  7.88  x  10  \ 

5  -  1.42;  the  lower  pair  is  for  v/yQ  -  1.57  x  10  ,  £  -  0.28. 


Normalised  growth  rate  for  the  negative  mass  instability 
versus  v/yq  for  the  "base  case"  parameters:  a  •  0.5m,  b  ■  2.2m, 

-  0.99m,  rj  ■  1.01m,  y  ■  3,  £  ■  1.  A  solid  curve  denotes  data 
obtained  from  a  numerical  solution  of  the  eigenvalue  problem;  a 
dashed  curve  denotes  data  from  the  dispersion  relation,  Eq.  (31). 

Normalized  growth  rate  for  the  negative  mass  instability  versus 
layer  thickness.  All  parameters  take  their  base  case  values,  (see 
text)  except  that  the  peak  density  in  the  parabolic  profile 
is  10^/cc. 

Normalized  growth  rate  for  the  negative  mass  Instability  versus 

outer  wall  position.  All  parameters  take  their  base  case  values. 

(See  text.)  Near  synchronous  conditions,  £u>  «. u>  where  u  is  a 

o  v  v 

waveguide  mode  satisfying  b+  +  b_  «=  0,  the  negative  mass  instability 
is  strongly  enhanced  and  Eq.  (31)  should  be  solved  as  a  cubic 
polynomial.  The  synchronous  or  enhanced  negative  mass  instability  is 
often  called  the  cyclotron  maser  Instability. 

Normalized  growth  rate  for  the  negative  mass  instability  versus 
inner  wall  position.  All  parameters  take  their  base  case  values. 

(See  text.) 

Normalized  growth  rate  for  the  negative  mass  instability 
versus  yq*  All  parameters  take  their  base  case  values,  (see  text) 
except  that  the  peak  density  in  the  parabolic  profile  is  !0®/cc. 


10.  Normalized  growth  rate  for  the  negative  mass  instability  versus 


azimuthal  mode  number  £.  All  parameters  take  their  base  case 


values.  (See  text.)  A  +  denotes  a  result  from  Eq.  (31)  and 


an  x  denotes  a  result  from  a  numerical  solution  of  the  eigenvalue 


problem. 


Fig.  B-l .  Normalized  admittances  vs.  outer  wall  position  for  various  inner 


wall  positions  for  l  ■  1  and  me  (yQ-l)  ■  1  keV. 


Fig.  B-2.  As  in  Fig.  B-l  but  for  me  (y  -1)  =  5  keV. 


Fig.  B-3.  As  in  Fig.  B-l  but  for  me  (y  -1)  *  40  keV. 


Fig.  B-4.  As  in  Fig.  B-l  but  for  me  (y  -1)  *  300  keV.  The  zeroes  of  b+  +  b_ 


correspond  to  waveguide  cutoff  frequencies. 
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FIGURE  B-3 
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APPENDIX  C 


Program  ASTRON 


Program  ASTRON 


Input  data  file:  FOR015.DAT  (free  format) 

All  quantities  in  MKS  units 


100 

200  a ,  b 
300  r^ ,  r2 

400  n00'  pwr 
500  0o(b),  Boo 

600  KBR 

700  ,  e2 

800  II 

900  IEV 

l000 

L100  ZW,  L 
L200  KT1,  KT2,  KT3 

L300  TOL,  e3 

L400  12,  13 


Label 

Inner,  outer  cylinder  radii 

Inner,  outer  beam  radii 

Peak  number  density,  power  (Note  1) 

Bias  potential  at  r  =  b;  Applied  B  field 

(+1,-1)  =  (greater,  lesser)  of  two  roots 
of  eqn.  of  motion 

Error  allowed  in  calculating  U(r^);  HQ(ri) 

If  =  1,  write  equilibrium  quantities 
to  file  FOR# 2,0 .  DAT 

(0,1)  =  (stop,  continue)  after  calculation 
of  equilibrium 

Label 

Initial  guess  for  u)/o)(r, ),  azimuthal 
mode  number 

Number  of  calls  to  integration  package 
in  (a-*r^)  ,  (r^-*-^)  ,  (r2-*-b)  respectively 

Integration  error  control;  error  control 
for  eigenvalue  convergence  test 

12  =  1  => Write  eigenfunction  to  file 
FOR# 2 2. DAT 

13  =  1  =>  Write  (w-iu)0  (r) ) /u)Q  (r^)  and  d  to 

file  FOR# 2 4 .DAT. 


Note  1:  Functional  form  of  density  is  taken  as 


_ 

n  =  n 

4 (r2-r) (r-r1) 

pwr 

oo 

(r9-r, ) 

Output  data  files:  (All  binary) 
Equilibrium  quantities: 


FOR020.DAT:  r,  EQ,  Ho#  u,  £,  h,  me '/ {u>0t ' 

Eigenfunction: 

FOR022.DAT:  r,  Re(rEQ1),  Im(rEel)  for  r  =  a 

Singularity  locations: 

FOR024.DAT:  r,  Re  (  <cu-  £ojq  (r) )  /uQ  (r.^ )  ,  Re(d) 


)  for  r 

to  b. 

,  lm(d) 


t  v i  i  ([  tr  one  on,.  1  or 

C  COMMON  CODE  FOR  FROGRAM  AST  RON 

IMPLICIT  REALte(A-H.O-R.T-Y) 

IMPLICIT  COMPLEX! 1 fc  2 
P AR AME  TER (Nl-200) 

C  N1  MUST  BE  EVEN 

COMMON/OEOM/A.B.RJ  .RI-.R2R1 
COMMON/CNTRL/N. DR. DR2.FBR.EPS) ,EP52, IEV 

common/den/deno, fwp. 122,000 

COMMON/BIAS/PHIOD,  BOO 

COMMON / CONS T / OE . E M , E PSO .  X MUO . C , 0EMC2 

COMMON /F  IElDS ,'EO  ( 0:  N1  >  ,H0(0‘.N1  ) 

COMMON /PARTS /U  <0:  NI  ) 

COMMON / DVDHM / DM  <  O : N 1 > 

COMMON /NORM/ BO  1 . B012, DRH01 . DRH02. DRH03.  DDRHO 

COMMON /CNTRL2 /) T 1 .FT2. )  T2, T0L.EPS2 

COMMON /EV/7L-E-. EL 7 

COMMON / ADM  n  S / 2  BF  LUS . 2  BM I NUS 

COMMON /COEF / 2MAT («,«> 

COMMON/EOUIL/XI ( 0 :  N 1  > ■ HG  < 0 : N 1 ) ,W0B(0: N1 ) 

COMMON/OUTS/1 1,  12.  13 
DATA  OE.EM.EPSO.XMUO.C 

1  /  1.6023ti-l 9.  R.  109SD-21 . 6.  8543D-  12.12.  Sfcfc4Ii-7, 2.  9979DB/ 

DATA  0EMC2/ 1  -  957 1  D-fc/ 

DATA  ZMAT/16* (O.DO.O.DO)/ 

< 


Vj  o  t»T  ,t»  to  f  1 


F • DGRAM  ASTRON 
INCLUDE  'ASTRONCDM.FOR' 

REUING  20 
N=N1 

CAUL  INFUT1 
CALL  ECAL 
CALL  HVCAL 
FAINT  JO 

10  FORMAT ( / ,  *  EQUILIBRIUM  CALCULATION  COMPLETED' , /> 

B01=U (0)  /DSGRT  <U  (0)  M2+1 .  DO) 

D012=B01*»2 

GOO=DSQRT  <U  <N/2) **2+1 .  DO) 

W0Rl=E<01*C/fil 
D0100  NN=0, N 
R=R1 +DFLOAT  (NN) *DR 
U2=U<NN)*»2 
GMA=DSORT (U2+1.D0) 

XI (NN) =GEMC2*EDEN(R) »R*GMA/U2 
HG(NN)=0EMC2*E0(Nn>*R/U2 
WOP(NN)=C*U(NN> / (GMA*F:»W0R1> 

IF C 11. EG. 1 ) THEN 

M=HG(NN) /gma 

GMA2=GMA*»2 

GMA4=GMA2»*2 

BETA2= 1 . DO- 1 . D0/GMA2 

DWDE=-  (BETA2+2.D0»H*X1  (NN)  /GMA4 )  / 

1  (GMA*BETA2* (1. DO- XI (NN) /GMA2+HG (NN) • *  2 > > 

WRITE  (20)  SNGL  (R)  .  SNDL  (EO(NN)  )  ,SNDL  (HO  (NN)  )  ,  SNGL  (U(NN)  )  , 

1  SNGL ( X 1 (NN) ) . SNGL  < H ) . SNGL ( DUDE  > 

END  IF 

100  CONTINUE 

XIAV=DSIMF (4) / (R2-R1 ) 

ENS=DSIMP (5) / (R2-R1 ) 

F  NORM*  DS I MR ( fc  > 

UAV=DSIMP(7) /FNORM 
U2AV=DSIMP  (e>  /FNORM 
WAV=DSIMP(9) /FNORM 
W2AV=DSIMF (10) /FNORM 
XNUC=0.5DO*BOI2*DSJMF(  J1)/(R1*»2) 

DELU=EiSORT  (12.  DO*  (U2AV-UAV*  »  2  >  >  /UAV 
DELW=DSORT (12. DO* (W2AV-WAV»*2> ) /WAV 
GAV= DSORT  ( 1 .  D0+UAVM2) 

FR INT20. XNUG , X I AV. ENS . UAV . GAV . DELU . DElW . WOR 1 

20  FORMAT C  NU/GAMMA=’ , 1FD10.3, /. 

1  AVERAGE  VALUE  OF  (PLASMA  FREO/CIRC  FREO ) *  * 2= * , OFF  1 0 . 4 . 
'  AVERAGE  VALUE  OF  SELF  FIELD  INDEX  <N£) - ' , F10. 4 , / , 

•  DENSITY  WEIGHTED  AVERAGE  OF  BETA*3AMMA= ' ,F10. 4, /. 

•  (CORRESPONDING  GAMMA= ', F 1 0. 4 ,')',/ , 

'  FRACTIONAL  LINEAR  MOMENTUM  SPREAO= • . F 1 0.  4 , / , 

'  FRACTIONAL  ANGULAR  FREDUENCY  SPREAD-  ', F 1 0. 4 ./ , 

'  CIRCULATION  FREQ  AT  R=R1= • . 1FD10. 3. *  RAD/SEC') 

CURRL=H0 (0) -HO (N) 

PRINT21 . CURRL 

21  FORMAT C  CURRENT  PER  UNIT  LENGTH- ' . 1 FT 1 C .  3 .  1  A/M'./> 
IFdEV.NE.  1 )  ODTO2O0 

CALL  INPUT 2 

CALL  SHOOTER 

Z  BPM=  Z  BPLUS + Z  BM 1 NUE 

PR I NT 30, ZBFLUS, ZBMINUS. ZBPM 

30  FORMAT (//,'  B*  ~ ' , 1P2D12. 4 . / . '  t-  =',2112.4./, 

1  '  5UM= ’,2012.4,/) 

IF (12. EG. 1 ) THEN 
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GMA=DSQRT  (U  (NN)  » t 2-t  1 .  DO) 

BET=U(NN> /GMA 
ZOM=ZW-EL*WOfc(NN) 

ZOMW=ZOM/WOE  (NN) 

RHO=  1 . DO+DFlOAT  <NN)  »DDRHO 
R*RHD»Rl 

ZT 1=  1 .  DO+HG  (NN)  »*2~ZDMU*»2 
ZT2=ZW»ZW»RH0*RH0»B012-El2 
ZT3=BET»B01»ZW*RH0~EL 
ZD=ZT  1  *  ZT2+X I  (NN)»2T3»«2 
S1=R 

S2=DREAL(ZOM> 

S3=DREAL(ZD) 

S4=DIMAG(ZD) 

WRITE (24)51 .S2. S3. £4 

CONTINUE 

ENDIF 

CONT 1  HUE 

STOP 

END 

SUBROUTINE  INPUT  1 
INCLUDE  ■astroncom.  fdf.  ■ 

READ (15,*) 

READ ( 1 5 .  *  >  A  >  B 
READ ( 15. » >Rl . RZ 
T22= ( (R2-R1 ) / 2 . DO ) »*2 
DR= (R2-R1 ) /DFLOAT(N) 

DR2=DR/2. DO 

READ  ( 1 5.  *  > DENO. FWR 

READ  < 1 5 . *  >  PH I OE . BOO 

READ < 15, * ) EBR 

READUS,  *)EFS1 ,EFS2 

READ(15.»)ll 

READ (15.*) IEV 

RETURN 

END 

SUEROUTINE  ECAL 
INCLUDE  '  ASTRONCOri.  FOR  ' 

ELBA=DLOD(E/A) 

Cl=(FHIOB*DSIMP(l>  >/ELBA 

E0(0)=-C1/R1 

F-O.DO 

R=R1 

DO 10  NN= 1 . N 

P=F  +  (EDEN (R) +EDEN (R+DR) *4. DO*EDEN (R+DRO ) *DRC / 
R=R+DR 

E0(NN)=-(C1+F ) /R 

RETURN 

END 

SUBROUTINE  HVCAL 
INCLUDE  ’ ASTRONCOM. FOR ' 

AAE=QE*EO (0) »R1 / (EM *C*C) 

B2A2=B»»2-A**2 

AAH1=OE»BOO»0.5DO* (R1+R2) / (EM*C) 

BETA1 =AAH1 /DSORT ( 1 , D0+AAH1 * » 2 ) 

HOOM=BOO/XMUO 

HOM=HOON+C» BETA 1»DS IMF  (3)  .'B2AC 
LH=0 

CONT INUE 
F  RINT2.  HOfi 

FORMAT ( '  HOM=  '.1FD10.3) 

AAh=OE*  XMJ0»H0N»R1 / (EM*C ) 


m  a  u  nu 


-  _  .  .  -  _  -  -.1 —  7* — .r  . — .  } — ,,  '  '■  '  ■  >,■  yi  _« 


CALL  UROOTSCUM, AAE.AAH) 

U(0)=UM 

C  INTEGRATION  OF  U  EQUATION 

0020  NN=0.N-1 
Y=U(NN) 

R=R1*DFL0AT  (NN)  *0R 
E=E0(NN) 

F1=UF'RM(R.  Y,E> 

E= (EO (NN) *E0 CNN* 1 > ) /2. 00 
F2=UF'RM (R+0R2,  Y*0R2»F  1  ,E> 

F2=UFRM(R*0R2. Y+0R2*F2.E> 

E=EO  CNN* 1 ) 

F4=UFRM(R*DR, Y+DR*F3,E> 

D  UlNN+1 >=Y*  (F1*2.D0»  (F2*F3)*F4)  *0R/6.D0 

CALCULATION  OF  U  ARRAY  COMPLETED 
CALCULATE  XI  AND  HD  ARRAYS 
00  22  NN-0 . N 
R=F:1  +DF LGAT  (NN)*DR 
U2=U(NN) *»2 
CMA=OSORT  (U2* 1.00) 

XI (NN) =0EHC2*EDEN(R) *R*GMA/Ul 
2  HG(NN)=0EMC2»E0CNN)»R/U2 

CALCULATE  SOLUTION  for  DVO/DHO- 
Gt1A=0S0Rl (U(0) **2*1.00) 

0H(0)=H0M»QE»XMU0»R1 / (EM*C< (DMA* *2) * (GMA+HD (0) ) > 

DOR I* DM (0) 

0024  NN=0,N-1 
00R0-00R1 
QQR1=Q0R (NN+1 ) 

24  0M(NN*1 )=OM(NN) »OEXF (0. 500* (DQR0+00R1 ) *0R) 

FS=H0M/H00M-0SIMP(2) / (H00M*B2A2> -1 . DO 
HDF5=H0M/H00M*C»DSIMf< 12) / (H00M*B2A2) 

C  ADJUST  VALUE  FOR  HO- 

HOMNE  W=HOM* ( 1 . OO-FS/HOFS) 

IF (OABS( (HOMNEW-HOM) /HOM) . LE . EPS2 ) GDT020 
LH=LH*1 

IF (LH.GT.20)ST0F  -HOM  LOOF  CONV  FAILURE1 

H0M=H0MNEW 

G0T01 

30  CONTINUE 

PRINT  150. LH 

15C  FORMATC  HVCA; _ FINAL  CONV  IN  , 12,'  ITERATIONS',/) 

HO ( 0) =HOMNEW 
R=Rl 

0040  NN=1,N 

VO=C*  U (NN- 1 ) /OSORT (U(NN-1 ) **2*1.00) 

V01=C»U(NN) /OSQRT (U (NN) * *2* 1 . 00) 

V02= (V0+V01 ) / 2 .  OC 
HO ( NN) =H0(NN-1 )*EFSO* 

1  (EOEN(R) »VO/R*EOEN(R*OR) * V01 / (R+OR) +4. 00»EOEN (R*0R2) 

2  V02/ (R+DR2) ) »0R2/3. 00 
R=R*  OR 

40  CONTINUE 

RETURN 
END 

FUNCTION  OSIMF (J) 

I NCL  UDE  ’ AS  TRONCOn .FOR’ 

0005=0.00 
EVENS=0.00 
00  101  NN=1,N-1,2 
R=R1 +DFLOAT (NN) *0R 
101  OOOS=ODOS+F (J,R> 

00102  NN=2, N-2. 2 
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R=R1 +DFLOAT  (NN)  tDR 
EVENS=EVENS*F  (J.R) 

DSIMF=  (F  (J.Rl >  *F  (J,R2>  »4. DOtODDS+2.  DO t EVENS)  *  DR/ 3.  DO 

RETURN 

END 

FUNCTION  F  (J.R) 

INCLUDE  '  ASTRONCOr,.FQR' 

NN= J I DNNT ( (R-Rl > /DR> 

I F  ( J.  GE .  6 ) OOT 0599 
GOTOOOO. 200. 300. 400.500)  ,  J 
CONT 1NUE 

F=EDEN  (F;>  tDLOO  (fi/trr 

C0T01 

CONTINUE 

VO=C*U(NN) /DSORT (U(NN) *»2+l .DO) 

F=EDEN(R>  »  (R*R-R*fc)  *VO*EF'SO/R 

GOT  01 

CONTINUE 

F=EDEN(F:>  *  (R*R-D»fc)  IEFSO/R 
G0T01 
CONTINUE 
F  =  XI (NN) 

G0T01 

CONTINUE 

F=XI  (NN)  7(1.  DO+U  (NN)  M2) 

G0T01 

CONTINUE 

IF  (J.EO.  1  DGOTOl  100 
IF ( J. ED. 12) GOTO 1 200 
F=0. DO 

IF  (  (R.LE.R1 )  .OR.  (F..GE.R2)  JGOTOl 

FF=R* ( (R2-R) * (R-Rl > ) MFTJR 

GOTO ( 600 . 700 . 800 . 900 . 1 000 > ,  J-E 

CONTINUE 

F=FF 

GOTO! 

CONTINUE 

F=FF»U(NN> 

GOT  01 

CONTINUE 

F=FF» (U (NN) »»t> 

GOTO] 

CONTINUE 

F=FF*WOfc(NN) 

G0T01 

CONTINUE 

F=FF» (WOB(NN) **2) 

G0T01 

CONTINUE 

F=R*XI  (NN)  *  (WOE  (NN)  M2) 

GOTO) 

CONTINUE 

F  = (R*R-B*D) * (-EFSO»EDEN(R> /R) »DM(NN) 

RETURN 

END 

FUNCTION  EDEN (R  > 

INCLUDE  ’  ASTRONCOF. . FOR  ’ 

CALCULATES  -R*RH0D/EF&0 
EDEN=0. DO 

IF  (  (R.LE.Rl)  .OR.  (F..GE.R2)  )GGT01 
FR- (R2-R) » (R-Ri ) /T22 
DENS=DENO* (Fh»*FWR) 

EDEN=Q£»R»DENS/EFSO 


no  cn  ii  m  r*  n 


1  RETURN 

ENTRY  X1DL0GN(R,GAM,W0B1> 

C  CALCULATES  X1*DL0G<DEN)/DL0G(R> 

XIDLOGN=0.  DO 

IF (  (R.LE.R1  >  .OR.  (R.GE.R2) >G0T02 
FR=  (R2-R)  *  (R-Rl )  /T22 
F1*FR**(PWR-1.D0) 

F2=R» (R1+R2-2. DO*R> /T22 
F3=  <WOB (N/2> /WOB1 >  *»2 
F4=G00/GAM 

XIDLOGN=Xl (N/2) *PUR*F 1*F2*F3*F4 

2  RETURN 
END 

SUBROUTINE  UROOTS (UM, AAE , AAH) 

C  CALCULATES  ROOTS  OF  RADIAL  FORCE  BALANCE  EON 

IMPLICIT  REAL >6  ( A-H. O- Y  ) 

IMFLICIT  COMPLEX* 16  2 
DIMENSION  A(S> »R(4> , Z (4) 

COMMON/CNTRl/N.  DR.  DR2 . K BR, EPS1 . EPS2, IEv 
PRINT  1, AAE. A AH 

1  FORMATC  AAE=  ' . 1PD10.3./. *  AAH=  ' ,  IF'Dl 0. 3,  / > 

A ( 1 ) = 1 . DO 

A(2)=-2.D0»AAH 

A ( 2 ) = AAH* »  2-AAE »  »  2 

A <41=0. DO 

A(5)=-AAE»*2 

CALL  ZRFOLY ( A. 4 . 2 . 1 ER) 

F'RINTlO.Z 

10  FORMAT (4<6X,  1FD10.3.3X. 1PD10.3. /> > 

C  TEST  ROOTS 

NACCEPT =0 
D0100  Ml, 4 

IF (DABS (DIMAG (211  1  ) > . GE. EPS1 ) GOTOlOO 
U=DREAL (2 (M ) 

E=DAES <  U*  U-AAE  »  DSORT (U»U+1 . DO) -AAH*  U  > 

IF  (E.GE.EPS1 >G0  T0100 
NACCEF  T=NACCEPT  + 1 
PR1NT20.U 

20  FORMAT (6X.1FD1 0.3) 

R (NACCEPT >=U 
100  CONTINUE 

GOTO (200. 20 1.202.  203. 203> .  ( NACCEF' T ♦ 1 ) 

00  FRINT30 

0  FORMAT (//,  '  NO  EQUILIBRIUM  STATE  EXISTS',//) 

STOP 

01  PR1NT31 

1  FORMmT(//'  WARNING. . .only  ONE  EQUILIBRIUM  STATE  EXISTS',//) 

UM=R  ( 1 > 

G0T0300 

02  IF  (KBR. EQ. 1 ) UM=DMAX 1 (R(1>,R<2>) 

IF (KBR. EQ. - 1 ) UM=DMIN1 (R(l) ,R(2) > 

G0T0300 

03  PRINT33. NACCEPT 

3  FORMAT (//,'  ERROR. . . ' , 1 1 , '  APPARENT  EQUILIBRIUM  STATES',//) 

STOP 

300  CONTINUE 

GAM=D5DRT  <UM»*2+1 .DO) 

E)  1N=5. 1  IDS* (GAM- 1. DO) 

PRINT 40,  UM.  DAM. E).  IN 

40  FORMATC  AT  R=R1S',/,'  U=  ' , 1FD10. 3. 3X , ' GAMMA= 

1  0FFB.3.3X, 'EF  IN=  '.1PD10.3.'  VOLTS',/) 

RETURN 

END 
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FUNCTION  UPRM(R.UU.E) 

INCLUDE  ' ASTRONCOM. FOR ' 

GMA=DSORT <UU*UU+1 .  DO) 

GMA2=GMA»*2 

WP2=OE*ED£N(R)  /  (EM»GMA»R> 

WO=C*UU/ (GMA*R> 

W02=W0»#2 

HH=QEMC2#E*R/ (UU#UU*GMA> 

UPRM=  (UU/R) *  (1 .  DO-  (WP2/  (GMA2IW02)  >-GMA2*HH)  /  (1 

RETURN 

END 

REAL*6  FUNCTION  OOR(NN) 

INCLUDE  * ASTRONCOM . F  OR ' 

GMA=DSORT (U(NN) *»2+ 1 .DO) 

GMA2=C-MA*»2 
GMA4=GNA2»»2 
BETA2= 1 . DO-1 .  D0/GMA2 
H=HG(NN) /GMA 
H1=H* 1 . DO 
X  =  X1 (NN) 

R=Rl+DFLOAT (NN) *DR 

QOR=  <  l.D0-3.D0»BETA2>  *  ( 1 . DO- X/GMA2-GMA2»H)  /HI 

1  +  (X* (2. D0/GMA2- 1 . DO- A . DO* PETA2*H ) 

2  +HI (GMA4+4.D0*GMA2-1.D0-2.D0#BETA2»GMA4*H> ) 

2  / (GMA2» <H1*#2)  ) 

OQR=OOR/R 

RETURN 

END 

SUBROUTINE  1NRUT2 
INCLUDE  'ASTRONCOM. FOR' 

READ < 15.#) 

READ(15.*)ZW.L 

EL=DFLOAT (L) 

EL2=EL#*2 

READU5.  #)>  T1.KT2.)  T3 
DRH01= (Rl-A) / (DFLOAT (KT1 ) tRl ) 

DRH02= (R2-R1 ) / (DFLOAT (I  T2>  »R1 > 

DRH02= (B-R2) / (DFLOAT  O T3>  *Rl ) 

DDRHO=DR/Rl 

R2R1=R2/R1 

READ ( 15»  #  > TOL . EFS3 

READ(15.»)12.I3 

RETURN 

END 

subroutine  shooter 

IMPLICIT  REAL*8<A-H,0-R.T-Y> 

IMPLICIT  COMP  LEX  * lfc  Z 
COMMON/ GEOK/ A . B . R 1 . R2 . R2R 1 

COMMON /NORM /BO 1  ,  BO 1 2 , DRHOl . DRH02 . DRH03 . DDRHO 

C0MM0N/CNTRL2/E  T1 , RT2. K T3 , TOL . EPS 3 

COMMON /EV/ZW, EL. EL2 

COMMON / ADM ITS/ZBFLUS.ZRMINUS 

DIMENSION  V(B).WW(8,9> 

EXTERNAL  vf-rm, STORE 
PRINTS, JIDNNT(EL) 

5  FORMAT ( // '  EV  CALCULATION:  L='.I3,/> 

LL=0 

LZW=0 

1  PRINTIOC.  Zk'  LZW 

100  FORMAT <  ZW=  '.IF2DI2.4.'  ITERATION  ’.12) 

IF(LL.EO.  DG0T0200 

IF  (LZL'.CT.  15.' STOP  1  EV  LOOP  CONVERGENCE  FAILURE 
D010  1=1.6 


DO+HH) 
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VO.)  =0.I<0 
v  (3) =1 .  r»o 
RH0=A/R1 
REWIND  22 
CALL  STORE (V.RHD) 

CALL  STEF'PER  <8 .  V.  RHD>  DRH01 , FT  1 , TOL  >  WW ,  VPRM, STORE .  100) 
ZPMINUS=-EL*DCMPLX  <V(3)  .V(4) )  /DCMPLX  ( V  ( 1 )  .  V  (2)  ) 

RHO=  1 . DO 

CALL  STEPPER ( 6 . V . RHO , DRH02 . ) T2 . T  OL . WW , VPRM , STORE . 50 ) 
ZPPLUS=£L*DCMPLX(V(3> ,V(4> > /DCMPLX  < V ( 1 ) , V<2> ) 

RH0=R2Rl 

CALL  STEPPER! 8, V. RHO, DRH03,) T3. TOL . WW , VFRM, STORE , 100) 
ZUNEW=ZW-DCMPLX  (V(  1 )  ,  V(2>  >  /PCMFLX  <V<5>  .  V(6)  ) 

EREAL=DREAL  (  ZWNEW-ZW) /UREAL (ZW) 

EIMAG^O.DO 

IF  (DAPS  (DIMAS  <ZW)  )  .OE.  1 .  P-£>  E IMAGrPIMAG  (ZWNEW-ZW )  /PIMAG(ZW) 

IF  (DMAX  1  (DAPS  <EREAL)  .  DAPS  (El MAG )  )  .  LE. EPS3)  LL  =  1 

LZW=LZW<  1 

ZW=ZWNEW 

GDT01 

CONTINUE 

RETURN 

END 

SUPROUT INE  VPRM  < JX . RHO, V, DV) 

IMPLICIT  REAL >8 ( A-H, 0-R , T-Y > 

IMPLICIT  COMPLEX* lfc  Z 
COMMON/COEF/ZMAT (4.4) 

DIMENSION  V ( 6 ) • D V ( 8 ) 

DIMENSION  ZV <4 ) , ZDV ( 4 ) 

D010  J=l,4 

ZV(J)=DCMPLX  <V(2*U-1 ) . V(2* J>  > 

CALL  MATCAL (RHO) 

D020  J=l,4 
ZDV ( J) =0. DO 
D030  K-1,4 

ZDV(J)  =  ZDV(J)*ZMAT(J,K)»ZVO  ) 

CONTINUE 
D040  J=l,4 
J1=2»J-1 
J2=2*J 

DV( J1 ) =DREAL (ZDV ( J> > 

DV(J2)=DIMAG(ZDV(J>  > 

RETURN 

END 

SUPROUT INE  MATCAL (RHO)  * 

INCLUDE  ' ASTRONCOM. FOR ’ 

Y(Y1, Y2, Xl.DX,  X)=Y1+ (Y2-Y1 ) KX-X1 ) /PI 
IF ( (RHO.LE. l.DO) .OR. (RHO. GE . R2R1 ) ) GOTO 100 
NN1= JIOINT ( (RHO- l.DO) /DDRHO) 

NN2=NN1*1 

RHO 1 =  NN 1 »  DDRHO+ 1 . DO 

XI I=Y(X1 (NN1 )  ,  XI (NN2) ,RH01 , DDRHO, RHO) 

WOPI  =  Y (WOP(NNl)  .  WOP ( NN2 ), RHO 1, DDRHO, RHO) 

UI=Y (U (NN1 ) • U ( NN2 ) , RHO 1 . DDRHO . RHO) 

GAM2=1.D0+UI»*2 
GAM=DSORT (GAM2) 

PET^UI/GAM 

HG I =Y (HG (NN1 ) , HG (NN2 ) ,RH01 , DDRHO, RHO) 

H=HGI /GAM 
F  G2= 1 . DO+H 
ZOM=  ZW-EL*WOP  2 
Z  OMW=  Z  DM / WOP I 
Z  T 1  =  1 .  DO+HG I  *  *  2- Z OMW  *  1 2 
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ZWRC=ZW*B01 *RHO 

ZWRC2=ZWRC**2 

ZT2=ZWRC2-EL2 

ZT3=BET»ZWRC-£L 

ZDD=ZT1*ZT2+XI1*ZT3**2 

ZA=ZT1/ZDD 

T  =  <  1  .DO-X  1 1  /GAM2-H*GAH2> /  ( 1 .  DO+H) 

V0F-T/GAM2 
WOP=VOP- 1 .  DO 
GAMF-T»BET**2 

HP=-X1  1  /DAM2-H*  <3. D0*DAM2- 1 . DO)  *VOF 
XIEf-XIPLOGN  (RHDtRl  .  GAM, WOB) ) 

X1F-X1D-X1 1* ( OAMF  +2. DO* WOF) 

ZDDF-2. DO*  <GAM2»H» <OAMP*H+HF  >  +WDP*ZW*ZQMW/WOBI )  *ZT2 
♦2.D0*ZT1*ZWRC2+XIP*ZT3*»2 
+2.D0*XII*2T3*ZWRC»BET»  (l.DO+UOF ) 

ZDEFi=  X. I D+  X 1 1 »  <  -UOF  -OAMP+HF/  FG2-  Z  DDF  /  Z  DD ) 

ZC1= (ZOM/ZW-1 .  D0/GAM2)  4ZDER+2.D0*  *n»BET**2 
ZC=ZW*FG2»ZC1/  <ZOM»ZDD> 

♦XII* <ZT2-X1I»BET**2> / (ZDD*GAM2) 

♦  1 .  DO 

ZDDDW=-Z0MW*ZT2/W0B1+ZT1*ZWRC2/ZW 

♦  XII *BET*B01 *RHO*  2T3 
Z  DDDW=  2 . DO  *  ZDDDW / Z  DD 

ZDAW=-2 . DO*  ZOMW / (WOB I  *  Z DD ) -Z A*  Z BDDW 
Z0DF'W=2.D0»W0F»ZT2* (ZW+ZOM) / <WOBI»*2> 

♦4. DO*  <GAM2*H* <GAMP*H+HF> +WOP* ZW* ZOMW/WOBI ) *ZWRC2/ZW 
♦4 .  DO*  ZWRC2*  ZT 1  /ZW-4. DO*  ZWRC2* ZOMW/WOBI 
♦2.D0*XIF*ZT3»BET*B01*RHD 

♦2 . DO*  X 1 1 »  BET  »  BO 1 *RHO» ( VOP+  I . DO ) » ( 2 . DO*  BET  *  BO 1 »  Z W»  RHO-EL ) 
ZDFDW=  <ZDDPU-ZDC)DW*ZDDP>  /ZDD 
Z  DCW= -EL  *  WOB I  *  PG2*  ZC 1 / ( Z  OM»  ZOM*  Z  DD > 

-ZW*PG2*ZPDDW*2CI / (ZOMtZDD) 

♦ (ZW/ZOM) * (PG2/ZDD) * (ELtWOBI »ZDER/ (ZW**2) 

-XII* (ZOM/ZW- I.D0/GAM2) *ZDPDW) 

-XI I  * ZDDDW* ( ZT2-X l I*BET**2) / (ZDD4GAM2) 

♦  XII *2. DO*  ZWRC2/ <ZDD*DAM2*ZW> 

G0TQ200 

CONTINUE 
VACUUM  COEFS 

ZDEL=RH0*RH0»B012*ZW*ZW-EL2 
Z A= 1 . DO/ZDEL 

ZDAW=-ZA*ZA*2.D0*RH0*RH0*B012»ZW 
ZC=1 .DO 
ZDCW=0. DO 
CONTINUE 

ZMAT  < 1 , 2)  =  1 .  DO/ <  ZAJRHO) 

ZMAT (2. 1 >=-ZC/RHO 

ZMAT (3.2) =-ZDAW/ <  ZA*  Z A*RHD) 

ZMAT(3.4)=ZMAT<1 ,2) 

ZMAT<4,l)=-ZDCW/RHO 
ZMAT (4.3) =  ZMAT  <2- 1) 

RETURN 

END 

SUBROUTINE  STORE (V.RHO) 

IMPLICIT  REAL»B (A-H.0-R.T-Y) 

COMMON /GEOM/ A, B>  Rl . R2. R2R1 
COMMON 'OUTS/ 1 1 .12.13 
DIMENSION  V <B> 

IF  (12. NE .  DRETURN 
S','  1  =V  (  1  i 
S‘.'2=V  <2> 

SF=RHO*RJ 
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